实验四：求微分方程的解

常微分方程的数值解法专业班级：信息软件 姓名：吴中原 学号：120108010002 一、实验目的1、熟悉各种初值问题的算法，编出算法程序；2、明确各种算法的精度与所选步长有密切关系；通过计算更加了解各种 算法的优越性。

二、实验题目1、根据初值问题数值算法，分别选择二个初值问题编程计算；2、试分别取不同步长，考察某节点j x处数值解的误差变化情况； 3、试用不同算法求解某初值问题，结果有何异常； 4、分析各个算法的优缺点。

三、实验原理与理论基础（一） 欧拉法算法设计对常微分方程初始问题(6-1)(6-2)用数值方法求解时，我们总是认为(6-1)、(6-2)的解存在且唯一。

欧拉法是解初值问题的最简单的数值方法。

从(6-2)式由于y (x 0) = y 0已给定，因而可以算出),()('000y x f x y =。

设x 1 = h 充分小，则近似地有：),()(')()(00001y x f x y hx y x y =≈-(6-3)记 ,n ,,i x y y i i 10 )(== 从而我们可以取),(0001y x hf y y ==作为)(1x y 的近似值。

利用1y 及f (x 1, y 1)又可以算出)(2x y 的近似值：),(1112y x hf y y +=一般地，在任意点()h n x n 11+=+处)(x y 的近似值由下式给出),(1n n n n y x hf y y +=+(6-4)这就是欧拉法的计算公式，h 称为步长。

⎪⎩⎪⎨⎧==)( ),(d d 00y x y y x f x y（二）四阶龙格-库塔法算法设计：欧拉公式可以改写为：()111,i i i i y y k k hf x y +=+⎧⎪⎨=⎪⎩，它每一步计算(),f x y 的值一次，截断误差为()2o h 。

改进的欧拉公式可以改写为：()()()11212112,,i i i i i i y y k k k hf x y k hf x h y k +⎧=++⎪⎪=⎨⎪=++⎪⎩，它每一步要计算(),f x y 的值两次，截断误差为()3o h 。

常微分方程实验报告一、实验目的常微分方程是数学分析和实际应用中非常重要的一部分，本次实验的主要目的是通过实际操作和计算，深入理解常微分方程的概念、性质和求解方法，并能够将其应用到实际问题中，提高我们解决数学问题和实际应用问题的能力。

二、实验原理常微分方程是指含有一个自变量和一个未知函数及其导数的等式。

求解常微分方程的方法有很多，常见的有变量分离法、一阶线性方程的求解方法（如常数变易法）、恰当方程的求解方法（通过积分因子）等。

对于一阶常微分方程，形如＼（y' ＋ p(x)y ＝ q(x)＼）的方程，可以使用积分因子＼（e^｛＼int p(x)dx}＼）来求解。

对于可分离变量的方程，形如＼（g(y)dy ＝ f(x)dx\），可以通过分别积分求解。

三、实验内容（一）一阶常微分方程的求解1、求解方程＼（y' ＋ 2xy ＝ 2x\）首先，计算积分因子＼（e^｛＼int 2xdx} ＝ e^｛x^2}＼），然后将方程两边乘以积分因子得到：＼（（ye^｛x^2}）＇＝ 2xe^｛x^2}＼）两边积分可得＼（ye^｛x^2} ＝ e^｛x^2} ＋ C\），解得＼（y ＝1 ＋ Ce^｛x^2}＼）2、求解方程＼（xy' y ＝ x^2\）将方程化为＼（y' ＼frac{y}｛x} ＝ x\），这里＼（p(x) ＝＼frac{1}｛x}＼），积分因子为＼（e^｛＼int ＼frac{1}｛x}dx} ＝＼frac{1}｛x}＼）。

方程两边乘以积分因子得到＼（（＼frac{y}｛x}）＇＝ 1\），积分可得＼（＼frac{y}｛x} ＝ x ＋ C\），即＼（y ＝ x^2 ＋ Cx\）（二）二阶常微分方程的求解1、求解方程＼（y'＇ 2y' ＋ y ＝ 0\）特征方程为＼（r^2 2r ＋ 1 ＝ 0\），解得＼（r ＝ 1\）（二重根），所以通解为＼（y ＝（C_1 ＋ C_2x)e^x\）2、求解方程＼（y'＇＋ 4y ＝ 0\）特征方程为＼（r^2 ＋ 4 ＝ 0\），解得＼（r ＝＼pm 2i\），所以通解为＼（y ＝ C_1\cos(2x) ＋ C_2\sin(2x)＼）（三）应用常微分方程解决实际问题1、考虑一个物体在受到与速度成正比的阻力作用下的运动，其运动方程为＼（m\frac{dv}｛dt} ＝ kv\）（其中＼（m\）为物体质量，＼（k\）为阻力系数），求解速度＼（v\）随时间＼（t\）的变化。

《数学建模实验》实验报告学院名称数学与信息学院专业名称提交日期课程教师实验一：数学规划模型AMPL求解实验内容1. 用AMPL求解下列问题并作灵敏度分析:一奶制品加工厂用牛奶生产A1和A2两种奶制品，1桶牛奶可以在甲类设备上用12小时加工成3公斤A1或者在乙类设备上用8小时加工成4公斤A2，且都能全部售出，且每公斤A1获利24元，每公斤A2获利16元。

先加工厂每天能得到50桶牛奶的供应，每天工人总的劳动时间为480小时，并且甲类设备每天至多加工100公斤A1，乙类设备的加工能力没有限制，试为该厂制定一个计划，使每天的获利最大。

（1）建立模型文件：milk.modset Products ordered;param Time{i in Products }>0;param Quan{i in Products}>0;param Profit{i in Products}>0;var x{i in Products}>=0;maximize profit: sum{i in Products} Profit [i]* Quan [i]*x[i];subject to raw: sum{i in Products}x[i] <=50;subject to time:sum{i in Products}Time[i]*x[i]<=480;subject to capacity: Quan[first(Products)]*x[first(Products)]<=100;（2）建立数据文件milk.datset Products:=A1 A2;param Time:=A1 12 A2 8;param Quan:=A1 3 A2 4;param Profit:=A1 24 A2 16;(3) 建立批处理文件milk.runmodel milk.mod;data milk.dat;option solver cplex;solve;display x;(4)运行运行结果：CPLEX 11.0.0: optimal solution; objective 33602 dual simplex iterations (1 in phase I)x [*] :=A1 20A2 30;（5）灵敏度分析：model milk.mod;data milk.dat;option solver cplex;option cplex_options 'sensitivity';solve;display x;display x.rc, x.down, x.up;display raw, time, capacity;display raw.down, raw.up,raw.current, raw.slack;得到结果：【灵敏度分析】: x.rc x.down x.up:=A1 -3.55271e-15 64 96A2 0 48 72;raw = 48time = 2capacity = 0raw.down = 43.3333raw.up = 60raw.current = 50raw.slack = 0某公司有6个建筑工地，位置坐标为(a i, b i)(单位：公里),水泥日用量d i (单位：吨）1) 现有j j j吨，制定每天的供应计划，即从A, B两料场分别向各工地运送多少吨水泥，使总的吨公里数最小。

重庆大学
学生实验报告
实验课程名称
开课实验室_____________________________________ 学生姓名___________ 学号__________
开课时间2015 至2016学年第二学期
实验目的
总结与体会
设计记录表格，包括碰到的问题汇总及解决情况
遇到的问题解决方法
采用dsolvc函数求不出微分方程组的解该微分方程组较复杂，不能求出精确解。

可以采用Obe23
函数求出数值解
运行程序时出现错误调用obe23函数一定要注意格式
心得体会：采用MATLAB软件求解微分方程有多种方法，采用dsolve函数可以求山一些较简单微分方程的精确解，但是一些较复杂的微分方程或微分方程组无法求山其精确解，这吋我们可以采用obe23、obe45等函数求解微分方程（组）的数值解。

如果这种方法都无法求山，可以尝试向前欧拉法、向后欧拉法编程等。

注行距：选最小值16磅，每一图应有简短确切的题名，连IM）图号置于图下。

每一表应有简短确切的题名，连M表号置
于表上。

图农的题名及其中的文字采用小5号宋体。

公式应该有编号，编号靠右端。

教师签名
年月曰
备注：
1、同一章的实验作为一个实验项0，每个实验做完后提交电子稿到服务器的“全校任选课数学
实验作、Ik.提交”文件夹，文件名为“学院学号姓名实验几”，如“机械20073159张新实验一
2、提交的纸质稿要求双面打印，中途提交批改不需要封面，但最后一次需将该课程所有实验项
目内页与封面一起装订成册提交。

3、综合实验要求3人合作完成，请在实验报告上注明合作者的姓名。

第1篇一、实验背景随着科技的不断发展，数学实验在各个领域中的应用越来越广泛。

数学实验作为一种以计算机为工具，通过模拟、计算和验证等方法，对数学理论进行实践探索和研究的方法，已经成为数学研究的重要手段。

本次实验旨在通过数学实验，加深对数学理论的理解，提高数学应用能力，培养创新意识和团队协作精神。

二、实验目的1. 熟悉数学实验的基本方法，掌握数学实验的基本步骤。

2. 通过实验，加深对数学理论的理解，提高数学应用能力。

3. 培养创新意识和团队协作精神，提高自身综合素质。

三、实验内容本次实验主要包括以下内容：1. 实验一：线性方程组的求解通过编写程序，实现线性方程组的直接法、迭代法等求解方法，并对比分析各种方法的优缺点。

2. 实验二：矩阵运算实现矩阵的加法、减法、乘法、转置等基本运算，以及求逆矩阵、特征值和特征向量等高级运算。

3. 实验三：数值积分通过编写程序，实现定积分、变积分、高斯积分等数值积分方法，并分析各种方法的误差和适用范围。

4. 实验四：常微分方程的数值解法实现欧拉法、龙格-库塔法等常微分方程的数值解法，并对比分析各种方法的稳定性、精度和适用范围。

四、实验过程1. 确定实验内容，明确实验目的。

2. 设计实验方案，包括实验步骤、算法选择、数据准备等。

3. 编写实验程序，实现实验方案。

4. 运行实验程序，收集实验数据。

5. 分析实验数据，得出实验结论。

6. 撰写实验报告，总结实验过程和结果。

五、实验结果与分析1. 实验一：线性方程组的求解通过实验，验证了直接法和迭代法在求解线性方程组时的有效性。

直接法在求解大规模线性方程组时具有较好的性能，而迭代法在求解稀疏线性方程组时具有较好的性能。

2. 实验二：矩阵运算实验结果表明，矩阵运算的程序实现具有较高的精度和效率。

在实际应用中，可以根据具体需求选择合适的矩阵运算方法。

3. 实验三：数值积分通过实验，验证了各种数值积分方法的有效性。

高斯积分具有较高的精度，但在求解复杂函数时，需要调整积分区间和节点。

实验2--微分方程(基础实验)119 项目四 无穷级数与微分方程实验2 微分方程（基础实验）实验目的 理解常微分方程解的概念以及积分曲线和方向场的概念,掌握利用Mathematica 求微分方程及方程组解的常用命令和方法.基本命令1. 求微分方程的解的命令DSolve对于可以用积分方法求解的微分方程和微分方程组,可用Dsolve 命令来求其通解或特解.例如,求方程023=+'+''y y y 的通解, 输入DSolve[y ''[x]+3y '[x]+2y[x]==0,y[x],x]则输出含有两个任意常数C[1]和C[2]的通解:{}{}]2[C e ]1[C e ]x [y x x 2--+→注:在上述命令中,一阶导数符号 ' 是通过键盘上的单引号 ' 输入的,二阶导数符号 '' 要输入两个单引号,而不能输入一个双引号.又如,求解微分方程的初值问题:,10,6,03400='==+'+''==x x y y y y y输入Dsolve[{y''[x]+4 y'[x]+3y[x]==0,y[0]==6, y'[0]==10},y[x],x](*大括号把方程和初始条件放在一起*)则输出{}{}x 2x 3e 148(e ]x [y +-→-2. 求微分方程的数值解的命令NDSolve对于不可以用积分方法求解的微分方程初值问题,可以用NDSolve 命令来求其特解.例如要求方程5.0,032=+='=x y x y y的近似解)5.10(≤≤x , 输入NDSolve[{y'[x]==y[x]^2+x^3,y[0]==0.5},y[x],{x,0,1.5}](*命令中的{x,0,1.5}表示相应的区间*)则输出{{y->InterpolatingFunction[{{0.,1.5}},< >]}}注:因为NDSolve 命令得到的输出是解)(x y y =的近似值. 首先在区间[0,1.5]内插入一系 列点n x x x ,,,21Λ, 计算出在这些点上函数的近似值n y y y ,,,21Λ, 再通过插值方法得到 )(x y y =在区间上的近似解.3. 一阶微分方程的方向场一般地,我们可把一阶微分方程写为),(y x f y ='的形式,其中),(y x f 是已知函数. 上述微分方程表明:未知函数y 在点x 处的斜率等于函数120f 在点),(y x 处的函数值. 因此,可在Oxy 平面上的每一点, 作出过该点的以),(y x f 为斜率 的一条很短的直线(即是未知函数y 的切线). 这样得到的一个图形就是微分方程),(y x f y ='的方向场. 为了便于观察, 实际上只要在Oxy 平面上取适当多的点,作出在这些点的函数的 切线. 顺着斜率的走向画出符合初始条件的解,就可以得到方程),(y x f y ='的近似的积分曲 线.例如, 画出0)0(,12=-=y y dxdy 的方向场. 输入<<Graphics`PlotField`g1=PlotVectorField[{1,1-y^2},{x,-3,3},{y,-2,2}, Frame->True,ScaleFunction->(1&),ScaleFactor->0.16,HeadLength->0.01,PlotPoints->{20,25}];则输出方向场的图形（图2.1）, 从图中可以观察到, 当初始条件为2/10=y 时, 这个微分方程的解介于1-和1之间, 且当x 趋向于-∞或∞时, )(x y 分别趋向于1-与1.-3-2-10123-2-1012 -3-2-10123-2-112下面求解这个微分方程, 并在同一坐标系中画出方程的解与方向场的图解. 输入sol=DSolve[{y'[x]==1-y[x]^2,y[0]==0},y[x],x];g2=Plot[sol[[1,1,2]],{x,-3,3},PlotStyle->{Hue[0.1],Thickness[0.005]}];Show[g2,g1,Axes->None,Frame->True];则输出微分方程的解xxe e x y 2211)(++-=,以及解曲线与方向场的图形（图2.2）. 从图中可以看到, 微分方程的解与方向场的箭头方向相吻合.实验内容用Dsolve 命令求解微分方程例2.1 (教材 例2.1) 求微分方程 22x xe xy y -=+'的通解.输入Clear[x,y];DSolve[y '[x]+2x*y[x]==x*Exp[-x^2],y[x],x]或DSolve[D[y[x],x]+2x*y[x]==x*Exp[-x^2],y[x],x]则输出微分方程的通解:121 ⎭⎬⎫⎩⎨⎧⎭⎬⎫⎩⎨⎧+→--]1[C e x e 21]x [y 22x 2x 其中C[1]是任意常数.例2.2 (教材 例2.2) 求微分方程0=-+'x e y y x 在初始条件e y x 21==下的特解. 输入Clear[x,y];DSolve[{x*y ' [x]+y[x]-Exp[x]==0,y[1]==2 E},y[x],x]则输出所求特解:⎪⎭⎪⎬⎫⎪⎩⎪⎨⎧⎪⎭⎪⎬⎫⎪⎩⎪⎨⎧+→x e e ]x [y x 例2.3 (教材 例2.3) 求微分方程x e y y y x 2cos 52=+'-''的通解.输入DSolve[y ''[x]-2y '[x]+5y[x]==Exp[x]*Cos[2 x],y[x],x]//Simplify则输出所求通解:⎭⎬⎫⎩⎨⎧⎭⎬⎫⎩⎨⎧-++→])x 2[Sin ])1[c 4x (2]x 2[Cos ])2[c 81((e 81]x [y x 例2.4 (教材 例2.4) 求解微分方程x e x y +=''2, 并作出其积分曲线.输入g1=Table[Plot[E^x+x^3/3+c1+x*c2,{x,-5,5},DisplayFunction->Identity],{c1,-10,10,5},{c2,-5,5,5}];Show[g1,DisplayFunction->$DisplayFunction]; -4-224-40-20204060图2.3例2.5 (教材 例2.5) 求微分方程组⎪⎪⎩⎪⎪⎨⎧=--=++02y x dtdy e y x dt dx t 在初始条件0,100====t t y x 下的特解.输入122Clear[x,y,t];DSolve[{x' [t]+x[t]+2 y[t]==Exp[t], y'[t] -x[t]- y[t]==0,x[0]==1,y[0]==0},{x[t],y[t]},t]则输出所求特解:⎭⎬⎫⎩⎨⎧⎭⎬⎫⎩⎨⎧+-→→])t [Sin ]t [Cos e (21]t [y ],t [Cos ]t [x t例2.6 验证c y y x =+--)3305(15152是微分方程2)(42-='y x x y 的通解. 输入命令<<Graphics`PlotField`<<Graphics`ImplicitPlot`sol=(-5x^3-30y+3y^5)/15==C;g1=ImplicitPlot[sol/.Table[{C->n},{n,-3,3}],{x,-3,3}];g2=PlotVectorField[{1,x^2/(y^4-2)},{x,-3,3},{y,-3,3},Frame->True,ScaleFunction->(1&),ScaleFactor->0.16,HeadLength->0.01,PlotPoints->{20,25}];g=Show[g2,g1,Axes->None,Frame->True];Show[GraphicsArray[{g1,g2,g}]];则分别输出积分曲线如图 2.4(a), 微分方程的方向场如图 2.4(b). 以及在同一坐标系中画出积分曲线和方向场的图形如下图2.4 (c).-3-2-1123-2-112-3-2-10123-3-2-10123-3-2-10123-3-2-10123图2.4从图 2.4(c)中可以看出微分方程的积分曲线与方向场的箭头方向吻合, 且当∞→x 时, 无论初始条件是什么, 所有的解都趋向于一条直线方程.例2.7 (教材 例2.6) 求解微分方程,)1(122/5+=+-x x y dx dy 并作出积分曲线. 输入<<Graphics`PlotField`DSolve[y' [x]-2y[x]/(x+1)==(x+1)^(5/2),y[x],x]则输出所给积分方程的解为 ⎭⎬⎫⎩⎨⎧⎭⎬⎫⎩⎨⎧+++→]1[C )x 1()x 1(32]x [y 22/7123 下面在同一坐标系中作出这个微分方程的方向场和积分曲线(设),3,2,1,0,1,2,3---=C 输入t=Table[2(1+x)^(7/2)/3+(1+x)^2c,{c,-1,1}];g1=Plot[Evaluate[t],{x,-1,1},PlotRange->{{-1,1},{-2,2}},PlotStyle->RGBColor[1,0,0],DisplayFunction->Identity];g2=PlotVectorField[{1,-2y/(x+1)+(x+1)^(5/2)},{x,-0.999,1},{y,-4,4},Frame->True,ScaleFunction->(1&), ScaleFactor->0.16,HeadLength->0.01,PlotPoints->{20,25},DisplayFunction->Identity];Show[g1,g2,Axes->None,Frame->True,DisplayFunction->$DisplayFunction];则输出积分曲线的图形（图2.5）.-0.75-0.5-0.2500.250.50.751-1.5-1-0.50.511.52图2.5例2.8 求解微分方程,2)21(22-+='-y x y xy 并作出其积分曲线.输入命令<<Graphics`PlotField`DSolve[1-2*x*y[x]*y' [x]==x^2+(y[x])^2-2,y[x],x]则得到微分方程的解为.)2(323C y x x y ++-+= 我们在33≤≤-C 时作出积分曲线, 输入命令t1=Table[(3+Sqrt[3])Sqrt[3+24x^2-4x^4-4*c*x]/(6*x),{c,-3,3}];t2=Table[(3-Sqrt[3])Sqrt[3+24x^2-4x^4-4*c*x]/(6*x),{c,-3,3}];gg1=Plot[Evaluate[t1],{x,-3,3},PlotRange->{{-3,3},{-3,3}},PlotStyle->RGBColor[1,0,0],DisplayFunction->Identity];124gg2=Plot[Evaluate[t2],{x,-3,3},PlotRange->{{-3,3},{-3,3}},PlotStyle->RGBColor[1,0,0],DisplayFunction->Identity];g1=ContourPlot[y-x^3/3-x*(-2+y^2),{x,-3,3},{y,-3,3},PlotRange->{-3,3},Contours->7,ContourShading->False,PlotPoints->50,DisplayFunction->Identity];g2=PlotVectorField[{1,(x^2+y^2-2)/(1-2*x*y)},{x,-3,3},{y,-3,3},Frame->True,ScaleFunction->(1&),ScaleFactor->0.16,HeadLength->0.01,PlotPoints->{20,25},DisplayFunction->Identity];Show[g1,g2,Axes->None,Frame->True,DisplayFunction->$DisplayFunction];Show[gg1,gg2,g2,Axes->None,Frame->True,DisplayFunction->$DisplayFunction];则输出微分方程的向量场与积分曲线, 并输出等值线的图2.6.-3-2-10123-2-10123-2-10123-2-1123图2.6用NDSolve 命令求微积分方程的近似解例2.9 (教材 例2.7) 求初值问题:1,0)1()1(2.1=='-++=x y y xy y xy 在区间[1.2,4]上的近似解并作图.输入fl=NDSolve[{(1+x*y[x])*y[x]+(1-x*y[x])*y'[x]==0,y[1.2]==1},y,{x,1.2,4}]则输出为数值近似解(插值函数)的形式:{{y->InterpolatingFunction[{{1.2,4.}},< >]}}用Plot 命令可以把它的图形画出来.不过还需要先使用强制求值命令Evalu-ate, 输入 Plot[Evaluate[y[x]/.fl],{x,1.2,4}]则输出近似解的图形（图2.7）.125 1.5 2.53 3.5410203040图2.7如果要求区间[1.2,4]内某一点的函数的近似值, 例如8.1=x y ,只要输入y[1.8]/.fl则输出所求结果{3.8341}例2.10 (教材 例2.8) 求范德波尔(Van der Pel)方程5.0,0,0)1(002-='==+'-+''==x x y y y y y y在区间[0,20]上的近似解.输入 Clear[x,y];NDSolve[{y''[x]+(y[x]^2-1)*y'[x]+y[x]==0,y[0]==0,y'[0]==-0.5},y,{x,0,20}];Plot[Evaluate[y[x]/.%],{x,0,20}]可以观察到近似解的图形（图2.8）.5101520-2-112图2.8126 ⎪⎩⎪⎨⎧==+-'1)1(01sin 2y x y x y x 的数值解, 并作出数值解的图形.输入命令<<Graphics`PlotField`sol=NDSolve[{x*y'[x]-x^2*y[x]*Sin[x]+1==0,y[1]==1},y[x],{x,1,4}];f[x_]=Evaluate[y[x]/.sol];g1=Plot[f[x],{x,1,4},PlotRange->All,DisplayFunction->Identity];g2=PlotVectorField[{1,(x^2*y*Sin[x]-1)/x},{x,1,4},{y,-2,9},Frame->True,ScaleFunction->(1&),ScaleFactor->0.16,HeadLength->0.01,PlotPoints->{20,25},DisplayFunction->Identity];g=Show[g1,g2,Axes->None,Frame->True];Show[GraphicsArray[{g1,g}],DisplayFunction->$DisplayFunction];则输出所给微分方程的数值解及数值解的图2.9.1.522.533.544681 1.52 2.53 3.54-22468例2.11 (教材 例2.9) 求出初值问题⎪⎩⎪⎨⎧='==+'+''0)0(,1)0(cos sin 22y y xy x y y的数值解, 并作出数值解的图形.输入NDSolve[{y''[x]+Sin[x]^2*y'[x]+y[x]==Cos[x]^2,y[0]==1,y'[0]==0},y[x],{x,0,10}]127 Plot[Evaluate[y[x]/.%],{x,0,10}];则输出所求微分方程的数值解及数值解的图形（图2.10）.2468100.20.40.60.8图2.10例2.12 (教材 例2.10) 洛伦兹(Lorenz)方程组是由三个一阶微分方程组成的方程组.这三个方程看似简单, 也没有包含复杂的函数, 但它的解却很有趣和耐人寻味. 试求解洛伦兹方程组,0)0(,4)0(,12)0()(4)()()()()(45)()()()(16)(16)(⎪⎪⎩⎪⎪⎨⎧===-='-+-='-='z y x t z t y t x t z t y t x t z t x t y t x t y t x 并画出解曲线的图形.输入Clear[eq,x,y,z]eq=Sequence[x'[t]==16*y[t]-16*x[t],y'[t]==-x[t]*z[t]-y[t]+45x[t],z'[t]==x[t]*y[t]-4z[t]];sol1=NDSolve[{eq,x[0]==12,y[0]==4,z[0]==0},{x[t],y[t],z[t]},{t,0,16},MaxSteps->10000];g1=ParametricPlot3D[Evaluate[{x[t],y[t],z[t]}/.sol1],{t,0,16},PlotPoints->14400,Boxed->False,Axes->None];则输出所求数值解的图形（图2.11(a)）. 从图中可以看出洛伦兹微分方程组具有一个奇异吸引子, 这个吸引子紧紧地把解的图形“吸”在一起. 有趣的是, 无论把解的曲线画得多长, 这些曲线也不相交.128图2.11改变初值为,10)0(,10)0(,6)0(=-==z y x 输入sol2=NDSolve[{eq,x[0]==6,y[0]==-10,z[0]==10}, {x[t],y[t],z[t]},{t,0,24},MaxSteps->10000];g2=ParametricPlot3D[Evaluate[{x[t],y[t],z[t]}/.sol2],{t,0,24},PlotPoints->14400,Boxed->False,Axes->None];Show[GraphicsArray[{g1,g2}]];则输出所求数值解的图形（图2.11(b)）. 从图中可以看出奇异吸引子又出现了, 它把解“吸”在某个区域内, 使得所有的解好象是有规则地依某种模式缠绕.实验习题1. 求下列微分方程的通解:(1) ;0136=+'+''y y y(2) ();024=+''+y y y(3) ;2sin 52x e y y y x =+'-''(4) .)1(963x e x y y y +=+'-''2. 求下列微分方程的特解:(1) ;15,0,029400='==+'+''==x x y y y y y(2) .1,1,02sin ='==++''==ππx x y yx y y 3. 求微分方程0cos 2)1(2=-+'-x xy y x 在初始条件10==x y 下的特解.分别求精确解和数值解)10(≤≤x 并作图.4. 求微分方程组⎪⎪⎩⎪⎪⎨⎧=--=++t t e y x dt dy e y x dt dx 235的通解.129 5. 求微分方程组⎪⎪⎩⎪⎨⎧==+-==-+==4,081,0300t t y y x dt dyxy x dt dx 的特解. 6. 求欧拉方程组324x y y x y x =-'+''的通解.7. 求方程5,0,011='==+'+''==x x y y y y x y 在区间[0,4]上的近似解.。

微分方程求解一、实验目的与要求1．掌握用Matlab 求微分方程及其方程组解的方法；2．学会求微分方程近似解的欧拉折线法；3．学会建立一些简单问题的微分方程模型，并能运用Matlab 分析研究这些问题。

二、问题描述对于很多实际问题，要直接找出所需的函数关系往往非常困难，但根据实际问题所提供的条件，有时却可以列出含有未知函数导数的关系式，这样的关系式就是所谓的微分方程。

怎样利用微分方程求得所需未知函数，往往是我们解决实际问题经常需要面对的问题，即解微分方程。

这里我们借用Matlab 对此问题进行简单探讨。

三、问题分析在处理关于微分方程的实际问题时，我们一般须先建立微分方程，再利用所学的数学知识解微分方程。

事实上真正能找到精确解的微分方程只是很少一部分，大部分只能求近似解，即数值解。

四、试验过程1．求微分方程解析解的命令。

求微分方程解析解的命令为：dsolve（'方程1''方程2'…，’初始条件1''初始条件2'…，’自变量’）, 对于可用积分方法求解的微分方程和微分方程组，可以用dsolve命令来求其通解和特解。

例1：要求方程y 3y 4y 0的通解，可以输入以下语句Matlab 命令：dsolve ('D2y+3*Dy-4*y=0', 'x')ans =C1*exp(-4*x)+C2*exp(x)运行结果：即y C i eC 2e x4x推。

如果自变量没有选定，默认自变量为2 •求微分方程数值解。

求微分方程数值解命令为 ode45，ode23, ode15& 对于不可以用积分方法求解的微分方程初值问题，可以用ode45，ode23,ode15s 命令求特解。

例4:求微分方程y yx y x 2 , y x 0 3的近似解(0 x 4)可用下 面的命令：fun cti on f=odef un 1(x,y) f=y*x+y+xA2;注：求一阶用D 表示，二阶导数用D2表示，三阶导数用D3表示，以此类例2:解方程y4y 5y 0运行结果: dsolve ('D2y+4*Dy+5*y=0','x')ans =C1*exp(-2*x)*si n(x)+C2*exp(-2*x)*cos(x)即: y e 2x C 1 sin x C 2 cosx例3:解方程1x 2 y 2xy xe"运行结果:即:dsolve ('(1+x A 2)*Dy+2*x*y=x*exp(x A 2)','x')ans=(1/2*exp(xA2)+C1)/(1+xA2) 1 x 2 -e x C 1 y 2 — y 1 x 2如果要求微分方程的初值问题：y 4y 2y 0 y x 0 6 , y 输入以下语句dsolve ('D2y+4*Dy-2*y=0','y(0)=6','Dy(0)=10','x')运行结果： ans =x010，可(3+11/6A(1/2))*exp((-2+6A(1/2))*x)+(-11/6A(1/2)+3)*exp(-(2+6A(1/2))*x)即:y 3 11 w611 2 6 x3 — e6[x,y]=ode45('odefu n1',[0,4],3);plot ( x , y , ' r -')输出结果为：例5:解初值问题y y t 1 , y 0 1dsolve('Dy=-y+t+1','y(0)=1') 输出结果：ans =t+exp(-t)即: y t e t现在我们用数值求解命令求解后和解析解比较function f=odef un 2(t,y)f=-y+t+1;t=0:0.1:1; y=t+exp(-t); plot(t,y,'b-')hold on[t,y]=ode45('odefu n2',[0,1],1);plot(t,y,'r.')hold off输出结果:例6:求初值问题y y sin 2x 0 , y 1 , y 1解：设y i y , y2 y，则原方程可化为y i y2y2y1sin 2xy i 1 , y2 1 Matlab 语言：fun cti on f=odef un 3(x,y)f=[y(2)；-y(1)-si n(2*x)];[x,y]=ode45('odefu n3',[pi,2*pi],[1,1]);plot(x,y(:,1),'r-')输出结果:利用ode45命令还可以求解耦合微分方程，所谓耦合微分方程，方程组中的未知函数是相互影响的，相互依赖的，其中的一个求解会影响到另一个求解，下面求一对耦合微分方程的数值解：例7:解方程组xt yt其中x 0 0 , y 0 2.1y t 0.01y t sin x tMatlab 语句：function f=odef un 4(t,y)f=[y (2),-0.01*y (2)-si n( y(1))]';[t,y]=ode15s('odefu n4',[0,100],[0,2.1]);函数x x t的图像：plot(t,y(:,1),'r-')输出结果为:函数y y t的图像：plot(t,y(:,2),'r-') 输出结果：用xt , yt生成参数图形plot(y(:,1),y(:,2),'r-') 输出结果:2iQ 5 10 15 2D 253 •欧拉折线法对于初值问题y f x, y , y X。

本科生实验报告实验课程微分方程数值解学院名称管理科学学院专业名称信息与计算科学学生学生学号指导教师林红霞实验地点6C402实验成绩二〇一五年十月二〇一五年十一月填写说明1、适用于本科生所有的实验报告〔印制实验报告册除外〕；2、专业填写为专业全称，有专业方向的用小括号标明；3、格式要求：①用A4纸双面打印〔封面双面打印〕或在A4大小纸上用蓝黑色水笔书写。

②打印排版：正文用宋体小四号，1.5倍行距，页边距采取默认形式〔上下2.54cm，左右2.54cm，页眉，页脚〕。

字符间距为默认值〔缩放100%，间距：标准〕；页码用小五号字底端居中。

③具体要求：题目〔二号黑体居中〕；摘要〔“摘要”二字用小二号黑体居中，隔行书写摘要的文字部分，小4号宋体〕；关键词〔隔行顶格书写“关键词”三字，提炼3-5个关键词，用分号隔开，小4号黑体)；正文部分采用三级标题；第1章××(小二号黑体居中，段前0.5行)1.1 ×××××小三号黑体×××××〔〕小四号黑体〔〕参考文献〔〕，参考文献用五号宋体，参照《参考文献著录规则〔GB/T 7714－2005〕》。

实验一常微分方程初值解法〔一〕1 实验内容分别用Euler法、改良Euler法、Runge-kutta法求解初值问题{u ·=−2∗t∗u2 u(0)=12 实验数据与实验结果1〕实验结果欧拉法的值与真实的值改良欧拉法的值与真实的值4阶Runge-kutta法的值与真实的值欧拉法的值与真实的值改良欧拉法的值与真实的值4阶Runge-kutta法的值与真实的值从上图可以看出，Euler法计算出的结果与其他方法所得结果相对来说精度太低。

4 程序代码清单实验二求解二阶常微分初值问题〔二〕1 实验内容分别用差分法、Jacobi迭代法求解初值问题{−dudx2=sin x 0<x<1u(0)=u(1)=02 实验数据与实验结果〔可用文字描述或贴图的方式进行说明〕1〕实验结果3 程序代码清单〔可直接将可运行源代码粘贴在下面的方框中〕实验三用向前向后差分格式求解〔三〕1 实验内容分别用向前向后格式求解初值问题{ðuðt=ð2uðu2(0≼x≼1,0<t≼1) u(x,0)=e x (0<x<1)u(0,t)=e t,u(1,t)=e1+t (0≼t≼1)2 实验数据与实验结果〔可用文字描述或贴图的方式进行说明〕1〕实验结果向前差分结果：T =Columns 1 through 10Columns 11 through 20Columns 21 through 30Columns 31 through 40Columns 41 through 50Columns 51 through 60Columns 61 through 70Columns 71 through 80Columns 81 through 90Columns 91 through 100Columns 101 through 110Columns 111 through 120Columns 121 through 130Columns 131 through 140Columns 141 through 150Columns 151 through 160Columns 161 through 170Columns 171 through 180Columns 181 through 190Columns 191 through 200Column 201x =Columns 1 through 10Column 11向后差分结果：图3-1 向后差分结果图3-2 向后差分结果4 程序代码清单〔可直接将可运行源代码粘贴在下面的方框中〕实验四显格式求解初值解法〔四〕1 实验内容用显格式求解{utt=uxx (0<x<1,t>0) u(0,t)=u(1,t)=0 (t>0)u(x,0)=sin4πx,ut(x,0)=sin8πx (0<x<1)2 实验数据与实验结果〔可用文字描述或贴图的方式进行说明〕1）实验结果图4-1 (x,u(:,501),'b')图图4-2 (x,u(:,1001),'b')图图4-3 (x,u(:,1501),'b')图4 程序代码清单〔可直接将可运行源代码粘贴在下面的方框中〕。

引言概述:高等数学数学实验报告（二）旨在对高等数学的相关实验进行探究与研究。

本次实验报告共分为五个大点，每个大点讨论了不同的实验内容。

在每个大点下，我们进一步细分了五到九个小点，对实验过程、数据收集、数据分析等进行了详细描述。

通过本次实验，我们可以更好地理解高等数学的概念和应用。

正文内容:一、微分方程实验1.利用欧拉法求解微分方程a.介绍欧拉法的原理和步骤b.详细阐述欧拉法在实际问题中的应用c.给出具体的实例，展示欧拉法的计算步骤2.应用微分方程建立模型求解实际问题a.介绍微分方程模型的建立方法b.给出一个具体的实际问题，使用微分方程建立模型c.详细阐述模型求解步骤和结果分析3.使用MATLAB求解微分方程a.MATLAB求解微分方程的基本语法和函数b.给出一个具体的微分方程问题，在MATLAB中进行求解c.分析结果的准确性和稳定性二、级数实验1.了解级数的概念和性质a.简要介绍级数的定义和基本概念b.阐述级数收敛和发散的判别法c.讨论级数的性质和重要定理2.使用级数展开函数a.介绍级数展开函数的原理和步骤b.给出一个函数，使用级数展开进行近似计算c.分析级数近似计算的精确度和效果3.级数的收敛性与运算a.讨论级数收敛性的判别法b.介绍级数的运算性质和求和法则c.给出具体的例题，进行级数的运算和求和三、多元函数极值与最值实验1.多元函数的极值点求解a.介绍多元函数的极值点的定义和求解方法b.给出一个多元函数的实例，详细阐述求解过程c.分析极值点对应的函数值和意义2.多元函数的条件极值与最值a.讨论多元函数的条件极值的判定法b.给出一个具体的多元函数，求解其条件极值和最值c.分析条件极值和最值对应的函数值和意义3.利用MATLAB进行多元函数极值与最值的计算a.MATLAB求解多元函数极值与最值的基本语法和函数b.给出一个多元函数的具体问题，在MATLAB中进行求解c.分析结果的准确性和可行性四、曲线积分与曲面积分实验1.曲线积分的计算方法与应用a.介绍曲线积分的定义和计算方法b.给出一个具体的曲线积分问题，详细阐述计算过程c.分析曲线积分结果的几何意义2.曲线积分的应用举例a.讨论曲线积分在实际问题中的应用b.给出一个实际问题，使用曲线积分进行求解c.分析曲线积分结果的实际意义和应用价值3.曲面积分的计算方法与应用a.介绍曲面积分的定义和计算方法b.给出一个具体的曲面积分问题，详细阐述计算过程c.分析曲面积分结果的几何意义五、空间解析几何实验1.空间曲线的参数方程表示与性质a.介绍空间曲线的参数方程表示和性质b.给出一个具体的空间曲线，转化为参数方程表示c.分析参数方程对应的几何意义和性质2.平面与空间直线的位置关系a.讨论平面与空间直线的位置关系的判定方法b.给出一个具体的平面与空间直线的问题，判定其位置关系c.分析位置关系对应的几何意义和应用实例3.空间直线与平面的夹角和距离计算a.介绍空间直线与平面的夹角和距离的计算方法b.给出一个具体的空间直线和平面，计算其夹角和距离c.分析夹角和距离计算结果的几何意义总结:通过本次高等数学数学实验报告（二），我们深入了解了微分方程、级数、多元函数极值与最值、曲线积分、曲面积分以及空间解析几何的相关概念和应用。