常微分方程的求解与定性分析实验报告完整版

福建农林大学计算机与信息学院（数学类课程）课程实习报告课程名称：常微分方程课程实习实习题目：常微分方程数值求解问题的实习姓名：上官火玉系：信息与计算科学专业：信息与计算科学年级：2008学号：081152048指导教师：陈永雪职称：讲师2010 年12 月 1 日福建农林大学计算机与信息学院数学类课程实习报告结果评定目录1. 实习的目的和任务 (1)2. 实习要求 (1)3. 实习地点 (1)4. 主要仪器设备 (1)5. 实习内容 (1)5.1 用不同格式对同一个初值问题的数值求解及其分析 (1)5.1.1求精确解 (1)5.1.2用欧拉法求解 (3)5.1.3用改进欧拉法求解 (5)5.1.4用4级4阶龙格—库塔法求解 (7)5.1.5 问题讨论与分析 (10)5.2 一个算法不同不长求解同一个初值问题及其分析 (12)6. 结束语 (13)参考文献 (13)常微分方程课程实习1.实习的目的和任务目的：通过课程实习能够应用MATLAB软来计算微分方程（组）的数值解；了解常微分方程数值解。

任务：通过具体的问题，利用MATLAB软件来计算问题的结果，分析问题的结论。

2.实习要求能够从案例的自然语言描述中，抽象出其中的数学模型；能够熟练应用所学的数值解计算方法；能够熟练使用MATLAB软件；对常微分方程数值解有所认识，包括对不同算法有所认识和对步长有所认识。

3.实习地点数学实验室4.主要仪器设备计算机、Microsoft Windows XPMatlab 6.55.实习内容5.1 用欧拉方法，改进欧拉方法，4阶龙格—库塔方法分别求下面微分方程的初值dy/dx=y+x+1 y(0)=1 x∈[0，2]5.1.1求精确解首先可以求得其精确解为：y=3*exp(x)-x-25.1.1 程序代码：>> x=0:0.1:2;>> y=3*exp(x)-x-2>> plot(x,y,'b*-');>> Data=[x',y']y =Columns 1 through 91.0000 1.2155 1.4642 1.74962.0755 2.44622.86643.3413 3.8766Columns 10 through 184.47885.1548 5.91256.76047.70798.76569.9451 11.2591 12.7218Columns 19 through 2114.3489 16.1577 18.1672Data =0 1.00000.1000 1.21550.2000 1.46420.3000 1.74960.4000 2.07550.5000 2.44620.6000 2.86640.7000 3.34130.8000 3.87660.9000 4.47881.0000 5.15481.1000 5.91251.2000 6.76041.3000 7.70791.4000 8.76561.5000 9.94511.6000 11.25911.7000 12.72181.8000 14.34891.9000 16.15772.0000 18.1672>>00.20.40.60.81 1.2 1.4 1.6 1.825.1.2 用欧拉法求解程序如下：建立函数文件cwfa1.mfunction [x,y]=cwfa1(fun,x_span,y0,h)x=x_span(1):h:x_span(2);y(1)=y0;for n=1:length(x)-1y(n+1)=y(n)+h*feval(fun,x(n),y(n)); endx=x';y=y';在MATLAB输入以下程序：>> clear all>> fun=inline(' y+x+1');>> [x,y]=cwfa1(fun,[0,2],1,0.1);>> [x,y]>> plot(x,y,'r*-')结果及其图象：ans =0 1.00000.1000 1.20000.2000 1.43000.3000 1.69300.4000 1.99230.5000 2.33150.6000 2.71470.7000 3.14620.8000 3.63080.9000 4.17381.0000 4.78121.1000 5.45941.2000 6.21531.3000 7.05681.4000 7.99251.5000 9.03171.6000 10.18491.7000 11.46341.8000 12.87981.9000 14.44772.0000 16.182500.20.40.60.81 1.2 1.4 1.6 1.825.1.3用改进欧拉法求解：程序如下：建立函数文件cwfa2.mfunction [x,y]=cwfa2(fun,x_span,y0,h)x=x_span(1):h:x_span(2);y(1)=y0;for n=1:length(x)-1k1=feval(fun,x(n),y(n));y(n+1)=y(n)+h*k1;k2=feval(fun,x(n+1),y(n+1));y(n+1)=y(n)+h*(k1+k2)/2;endx=x';y=y';在MATLAB输入以下程序：>> clear all>> fun=inline(' y+x+1');>> [x,y]=cwfa2(fun,[0,2],1,0.1); >> [x,y]>> plot(x,y,'r+-')结果及其图象：ans =0 1.00000.1000 1.21500.2000 1.46310.3000 1.74770.4000 2.07270.5000 2.44230.6000 2.86130.7000 3.33470.8000 3.86840.9000 4.46851.0000 5.14221.1000 5.89721.2000 6.74191.3000 7.68581.4000 8.73931.5000 9.91391.6000 11.22241.7000 12.67871.8000 14.29851.9000 16.09882.0000 18.0987>>00.20.40.60.81 1.2 1.4 1.6 1.825.1.4 用4阶龙格—库塔求解程序如下：建立函数文件cwfa3.mfunction [x,y]=cwfa3(fun,x_span,y0,h)x=x_span(1):h:x_span(2);y(1)=y0;for n=1:length(x)-1k1=feval(fun,x(n),y(n));k2=feval(fun,x(n)+h/2,y(n)+h/2*k1);k3=feval(fun,x(n)+h/2,y(n)+h/2*k2);k4=feval(fun,x(n+1),y(n)+h*k3);y(n+1)=y(n)+h*(k1+2*k2+2*k3+k4)/6;endx=x';y=y';在MATLAB输入以下程序：>> clear all;>> fun=inline(' y+x+1');>> [x,y]=cwfa3(fun,[0,2],1,0.1); >> [x,y]>> plot(x,y, 'b+-')结果及其图象：ans =0 1.00000.1000 1.21550.2000 1.46420.3000 1.74960.4000 2.07550.5000 2.44620.6000 2.86640.7000 3.34130.8000 3.87660.9000 4.47881.0000 5.15481.1000 5.91251.2000 6.76031.3000 7.70791.4000 8.76561.5000 9.94511.6000 11.25911.7000 12.72181.8000 14.34891.9000 16.15772.0000 18.1671>>00.20.40.60.81 1.2 1.4 1.6 1.825.1.5 问题讨论与分析由以上数值分析结果绘制表格：x=0:0.1:2;y1 = [1.0000 1.2155 1.4642 1.7496 2.0755 2.4462 2.8664 3.3413 3.8766 4.4788 5.1548 5.9125 6.7604 7.7079 8.7656 9.9451 11.2591 12.7218 14.3489 16.1577 18.1672];>> y1=[1.0000 1.2155 1.4642 1.7496 2.0755 2.4462 2.8664 3.3413 3.8766 4.4788 5.1548 5.9125 6.7604 7.7079 8.7656 9.9451 11.2591 12.7218 14.3489 16.1577 18.1672];>> y2=[1.0000 1.2000 1.4300 1.6930 1.9923 2.3315 2.7147 3.1462 3.6308 4.1738 4.7812 5.4594 6.2153 7.0568 7.9925 9.0317 10.1849 11.4634 12.8798 14.4477 16.1825];>> y3=[1.0000 1.2150 1.2192 1.4631 1.7477 2.0727 2.4423 2.8613 3.8684 4.4685 5.1422 5.8972 6.7419 7.6858 8.7393 9.9139 11.2224 12.6787 14.2985 16.0988 18.0987];>> y4=[1.0000 1.2155 1.4642 1.7496 2.0755 2.4462 2.8664 3.3413 3.8766 4.4788 5.1548 5.9125 6.7603 7.7079 8.7656 9.9451 11.2591 12.7218 14.3489 16.1577 18.1671 ];>> plot(x,y1,'r+-')>> hold on,plot(x,y2,'b-')>> plot(x,y1,'r+-')>> hold on,plot(x,y3,'b-')>> plot(x,y1,'r+-')>> hold on,plot(x,y4,'b-')00.20.40.60.81 1.2 1.4 1.6 1.8200.20.40.60.81 1.2 1.4 1.6 1.8200.20.40.60.81 1.2 1.4 1.6 1.82由上表和图可以看出欧拉法误差最大，而改进欧拉和龙格—库塔方法误差相对较小，并且龙格—库塔方法误差最小且大部分值都跟精确值相同。

常微分方程的数值解法专业班级：信息软件 姓名：吴中原 学号：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 。

实验报告实验项目名称常微分方程的数值解法实验室数学实验室所属课程名称微分方程数值解实验类型上机实验实验日期2013年3月11日班级10信息与计算科学学号2010119421姓名叶达伟成绩实验概述：【实验目的及要求】运用不同的数值解法来求解具体问题，并通过具体实例来分析比较各种常微分方程的数值解法的精度，为以后求解一般的常微分方程起到借鉴意义。

【实验原理】各种常微分方程的数值解法的原理，包括Euler法，改进Euler法，梯形法，Runge-Kutta方法，线性多步方法等。

【实验环境】（使用的软硬件）Matlab软件实验内容：【实验方案设计】我们分别运用Euler法，改进Euler法，RK方法和Adams隐式方法对同一问题进行求解，将数值解和解析解画在同一图像中，比较数值解的精度大小，得出结论。

【实验过程】（实验步骤、记录、数据、分析）我们首先来回顾一下原题：对于给定初值问题：1. 求出其解析解并用Matlab画出其图形；2. 采用Euler法取步长为0.5和0.25数值求解（2.16）,并将结果画在同一幅图中，比较两者精度;3. 采用改进Euler法求解（2.16），步长取为0.5；4. 采用四级Runge-Kutta法求解（2.16），步长取为0.5；5. 采用Adams四阶隐格式计算（2.16），初值可由四级Runge-Kutta格式确定。

下面，我们分五个步骤来完成这个问题：步骤一，求出（2.16）式的解析解并用Matlab 画出其图形； ，用Matlab 做出函数在上的图像，见下图：00.51 1.52 2.53 3.54 4.550.511.522.533.5x 1015y=exp(1/3 t 3-1.2t)exact solution图一 初值问题的解析解的图像步骤二，采用Euler 法取步长为0.5和0.25数值求解(2.16),并将结果画在同一幅图中，比较两者精度;我们采用Euler 法取步长为0.5和0.25数值求解，并且将数值解与解析解在一个图中呈现，见下图：00.51 1.52 2.53 3.54 4.550.511.522.533.5x 1015Numerical solution of Euler and exact solutionexact solution h=0.5h=0.25图二 Euler 方法的计算结果与解析解的比较从图像中不难看出，采用Euler 法取步长为0.5和0.25数值求解的误差不尽相同，也就是两种方法的计算精度不同，不妨将两者的绝对误差作图，可以使两种方法的精度更加直观化，见下图：00.51 1.52 2.53 3.54 4.550.511.522.533.5x 1015Absolute error of numerical solution and exact solutionh=0.5h=0.25图三 不同步长的Euler 法的计算结果与解析解的绝对误差的比较 从图像中我们不难看出，步长为0.25的Euler 法比步长为0.5的Euler 法的精度更高。

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

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

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

对于一阶常微分方程，形如＼（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\）的变化。

常微分方程实验报告《常微分方程》综合性实验实验报告实验班级05应数（3）学生姓名江晓荣学生学号200530770314指导老师方平华南农业大学理学院应用数学系实验微分方程在数学建模中的应用及数值解的求法一、实验目的1．了解常微分方程的基本概念。

2．常微分方程的解了解析解和数值解。

3．学习、掌握MA TLAB 软件有关求解常微分方程的解析解和数值解的有关命令。

4. 掌握微分方程在数学建模中的应用。

二、实验内容1．用MA TLAB 函数dsolve 符号求解常微分方程的通解和特解。

2．用MA TLAB 软件数值求解常微分方程。

三、实验准备1．用MA TLAB 求常微分方程的解析解的命令用MA TLAB 函数dsolve 求常微分方程()(,,,,,,)0n F x y y y y y ''''''= （7.1）的通解的主要调用格式如下：S=dsolve('eqn', 'var')其中输入的量eqn 是改用符号方程表示的常微分方程(,,,2,)0F x y Dy D y Dny = ，导数用D 表示，2阶导数用D2表示，以此类推。

var 表示自变量，默认的自变量为t 。

输出量S 是常微分方程的解析通解。

如果给定常微分方程（7.1）的初始条件()00010(),(),,()n n y x a y x a y x a '=== ，则求方程（7.1）的特解的主要调用格式如下：S=dsolve('eqn', 'condition1 ',…'conditonn ','var')其中输入量eqn ，var 的含义如上，condition1,…conditonn 是初始条件。

输出量S 是常微分方程的特解。

2．常微分方程的数值解法除常系数线性微分方程可用特征根法求解、少数特殊方程可用初等积分法求解外，大部分微分方程无解析解，应用中主要依靠数值解法。

1. 分别用Euler 法和ode45解下列常微分方程并与解析解比较： (1) ,(0)1,13y x y y x '=+=<<Euler 法：function [t,y]=euler(Fun,tspan,y0,h) t=tspan(1):h:tspan(2); y(1)=y0;for i=1:length(t)-1y(i+1)=y(i)+h.*feval(Fun,t(i),y(i)); end t=t'; y=y';function f=Fun(x,y) % 常微分方程的右端函数 f=x+y;>> [x,y]=euler('Fun',[0,3],1,0.1)>> [x,y] ans =0 1.0000 0.1000 1.1000 0.2000 1.2200 0.3000 1.3620 0.4000 1.5282 0.5000 1.7210 0.6000 1.9431 0.7000 2.1974 0.8000 2.4872 0.9000 2.8159 1.0000 3.1875 1.1000 3.6062 1.2000 4.0769 1.3000 4.6045 1.4000 5.1950 1.5000 5.8545 1.6000 6.5899 1.7000 7.4089 1.8000 8.3198 1.9000 9.3318 2.0000 10.4550 2.1000 11.7005 2.2000 13.0805 2.3000 14.6086 2.4000 16.2995 2.5000 18.1694 2.6000 20.2364 2.7000 22.5200 2.8000 25.0420 2.9000 27.8262 3.0000 30.8988ode45：>> [x,y]=ode45('Fun',[0,3],1) ans =0 1.0000 0.0502 1.0528 0.1005 1.1109 0.1507 1.17460.2010 1.2442 0.2760 1.3596 0.3510 1.4899 0.4260 1.63610.5010 1.7996 0.5760 1.9817 0.6510 2.1838 0.7260 2.40740.8010 2.6544 0.8760 2.9264 0.9510 3.2254 1.0260 3.55351.1010 3.9131 1.1760 4.3065 1.2510 4.7364 1.3260 5.20561.4010 5.7172 1.4760 6.2744 1.5510 6.8810 1.6260 7.54061.7010 8.2574 1.7760 9.0359 1.8510 9.8808 1.9260 10.79742.0010 11.7912 2.0760 12.8683 2.1510 14.0351 2.2260 15.29862.3010 16.6664 2.3760 18.1466 2.4510 19.7478 2.5260 21.47962.6010 23.3522 2.6760 25.3764 2.7510 27.5641 2.8260 29.92812.9010 32.4820 2.9257 33.3694 2.9505 34.2796 2.9752 35.21343.0000 36.1711解析解：>> y=dsolve('Dy=x+y','y(0)=1','x') y =2*exp(x) - x - 1(2) 20.01()2sin(),(0)0,(0)1,05y y y t y y t ''''-+===<< Euler 法：function f=Fun(t,y)% 常微分方程的右端函数f=[y(2);0.01*y(2)^2-2*y(1)+sin(t)];>> [t,y]=euler('Fun',[0,5],[0,1],0.2) ode45：>> [t,y]=ode45('Fun',[0,5],[0,1])t =0 0.0001 0.0001 0.0002 0.0002 0.0005 0.0007 0.0010 0.0012 0.00250.0037 0.0050 0.0062 0.0125 0.0188 0.0251 0.0313 0.0627 0.0941 0.12550.1569 0.2819 0.4069 0.5319 0.6569 0.7819 0.9069 1.0319 1.1569 1.28191.4069 1.5319 1.6569 1.7819 1.90692.0319 2.1569 2.2819 2.4069 2.53192.6569 2.7819 2.90693.0319 3.1569 3.2819 3.4069 3.5319 3.6569 3.78193.90694.0319 4.1569 4.2819 4.4069 4.5319 4.6569 4.7427 4.8285 4.91425.0000y =0 1.0000 0.0001 1.0000 0.0001 1.0000 0.0002 1.0000 0.0002 1.00000.0005 1.0000 0.0007 1.0000 0.0010 1.0000 0.0012 1.0000 0.0025 1.00000.0037 1.0000 0.0050 1.0000 0.0062 1.0000 0.0125 1.0000 0.0188 1.00000.0251 0.9999 0.0313 0.9998 0.0627 0.9987 0.0941 0.9965 0.1253 0.99340.1564 0.9893 0.2786 0.9632 0.3966 0.9220 0.5085 0.8662 0.6126 0.79670.7072 0.7146 0.7908 0.6210 0.8620 0.5176 0.9198 0.4058 0.9632 0.28760.9915 0.1647 1.0043 0.0392 1.0013 -0.0869 0.9826 -0.2117 0.9485 -0.33310.8996 -0.4490 0.8365 -0.5578 0.7605 -0.6577 0.6725 -0.7471 0.5742 -0.8246 0.4669 -0.8889 0.3525 -0.9393 0.2327 -0.9748 0.1095 -0.9950 -0.0154 -0.9996 -0.1398 -0.9887 -0.2619 -0.9624 -0.3798 -0.9212 -0.4916 -0.8657 -0.5957 -0.7970 -0.6904 -0.7161 -0.7742 -0.6242 -0.8460 -0.5228 -0.9046 -0.4134 -0.9491 -0.2978 -0.9789 -0.1777 -0.9934 -0.0549 -0.9945 0.0300 -0.9883 0.1146 -0.9748 0.1985 -0.9543 0.28092. 求一通过原点的曲线，它在(,)x y 处的切线斜率等于22,0 1.57.x y x +<<若x 上限增为1.58，1.60会发生什么？function f=Fun(x,y) % 常微分方程的右端函数 f=2*x+y.^2;>> [x,y]=ode45('Fun',[0,1.57],0) x =0 0.0393 0.0785 0.1178 0.1570 0.1963 0.2355 0.2748 0.3140 0.3533 0.3925 0.4318 0.4710 0.5103 0.5495 0.5888 0.6280 0.6673 0.7065 0.7458 0.7850 0.8243 0.8635 0.9028 0.9420 0.9813 1.0205 1.0598 1.0990 1.1383 1.1775 1.2168 1.2560 1.2953 1.3345 1.3738 1.4130 1.4248 1.4367 1.4485 1.4604 1.4722 1.4840 1.4959 1.5077 1.5140 1.5203 1.5265 1.5328 1.5376 1.5424 1.5472 1.5519 1.5543 1.5567 1.5591 1.5614 1.5631 1.5647 1.5664 1.5681 1.5685 1.5690 1.5695 1.5700 y =0 0.0015 0.0062 0.0139 0.0247 0.0386 0.0556 0.0758 0.09920.1259 0.1559 0.1895 0.2266 0.2675 0.3124 0.3615 0.4152 0.4738 0.5378 0.6076 0.6841 0.7679 0.8601 0.9620 1.0751 1.2014 1.3434 1.5045 1.6892 1.9037 2.1557 2.4577 2.8282 3.3003 3.9056 4.7317 5.9549 6.4431 7.0116 7.6832 8.4902 9.4821 10.7170 12.3090 14.4551 15.9220 17.7080 19.9390 22.8164 25.6450 29.2282 33.9673 40.5910 44.9434 50.3088 57.1229 66.1087 74.3108 84.7123 98.4901 117.7875 124.9206 132.9699 142.1268 152.6415若x 上限增为1.58，1.60,则超出运算的范围，发生溢出。

《常微分方程》实验报告四专业信息与计算科学班级姓名学号实验地点实验室实验时间2015.12.11实验名称：机械振动中的单摆问题实验目的：1．与计算机软件Maple结合起来，学会求常微分方程组解析解的方法；2．会用Maple软件做简单的计算机仿真；3．能够理论联系实际，从理论结论中客观分析实践中的现象。

实验内容：（给出实验程序与运行结果）问题：如图所示为一单摆示意图，物块的质量为，摆线长为，如图建立坐标系，（1）建立系统的运动函数，即物块坐标随时间的函数、；（2）设=10kg,L =1m，系统作微摆动，存在空气阻力，空气阻力与速度成正比，比例系数为u，u=2，求在初始条件x[0]=0 , [0]=0 下的系统运动函数、；（3）当系统作微摆动，不存在空气阻力时系统的运动函数、；（4）合理运用Matlab分别画出（2）、（3）中确定的各个函数的图象；（5）对（2）、（3）所得结果进行分析。

根据图象说明物块是做阻尼运动、无阻尼运动，还是周期Array运动；并分析给定参量具体数值后对反映系统的运动规律有无影响。

分析：1．对内容（1）的分析：从物块所做运动，明显看到它做圆周运动，另外作受力分析，可以发现物体运动遵循牛顿第二运动定律。

2．对内容（2）的分析：与（1）相比，物块做圆周运动是一样的；不同的是，物块作微摆动时，受力分析后我们作一个近似，T Mg≈，从而对在x方向上平衡建立一个方程。

利用Maltab 软件求方程组的解析解：一般情况下，对高阶方程首先引进新参量对方程进行降阶，之后直接使用dsolve 命令直接求出方程组的解。

但是这里求出的解可能并不简单明了，故而必须求它的数值解。

3．对内容（3），与（2）类似，而且相对简单。

由于这里不再受空气阻力。

4.对内容（4），用Matlab 软件中简单的作图命令plot 分别对（2）、（3）中所得函数进行做图。

5.对内容（5）通过分析（2）、（3）所得函数的图象，判断物块分别做什么样的运动。

常微分方程的边值问题和本征值问题一、问题描述利用搜索法和弦割法，得到该常微分方程的本征值,再利用打靶法计算多个本征值。

二、解决方法（一)搜索法1。

先随便猜测k的一个试验值,程序中令k=12.由Numerov算法根据本题的条件，kn+1=kn=kn-1=k,s=0,得到yn+2,yn+1，yn间的迭代公式令con=（k＊h）^2/12yn+2=2*(1-5＊con）＊yn+1/（1+con）—yn3自己给定φ的初始条件，然后利用公式得到边界值φ(1)4。

然后以小的步长dk增加k值，这里令dk=1,每当φ（1)改变符号时,就将步长减半后倒退回来重复5.当步长小于所要求的容许误差时终止程序，此时的k值即为所求。

（二）弦割法1.随便猜测两个k值,这里令k0=1,k1=22.自己给定φ的初始条件，对两个k值分别利用上述公式进行迭代，得到边界值y1(1）和y2（1）。

3.比较y1（1)和y2(1)的绝对值大小。

若绝对值大,说明对应的k值距离本征值距离较远。

4.将（k0+k1）/2赋给k2，边界值绝对值小的对应的k值保持不变，边界值绝对值大的对应k值重新定位k2的值.5.重复进行实验，当y1(1)和y(2）的差的绝对值小于容许误差时终止程序.此时k1的值即为所求。

当搜索法和弦割法大致求出了一个本征值后，利用打靶法，调整k值再度进行搜索，得到多个本征值，绘出其中一个本征值对应的函数图像，观察其性质。

三、程序实现1.搜索法subroutine add(t,y0，y1） !利用子程序表示函数值的迭代implicit nonereal(8）：：t,h，con，y0,y1，y2integer:：i，nn=10000h=1。

0/ncon=(t＊h)＊*2/12do i=1,n-1y2=2＊（1—5*con）＊y1/（1+con）—y0 ！利用Numerov算法,得到迭代公式y0=y1 !向前迭代y1=y2end doreturnend subroutine addprogram zy3implicit nonereal(8)：：diffk，dk，yold,k，b0,b1integer::sb0=0.01 !取初始值，根据题目条件，令y0=y1,来保证x=0的位置导数为0 b1=0。

学生实验报告（4）一、实验综述. 归纳和学习求解常微分方程(组)的基本原理和方法；2. 掌握解析、数值解法，并学会用图形观察解的形态和进行解的定性分析；3. 熟悉MATLAB软件关于微分方程求解的各种命令；4. 通过范例学习建立微分方程方面的数学模型以及求解全过程；通过该实验的学习，使学生掌握微分方程(组)求解方法（解析法、欧拉法、梯度法、改进欧拉法等），对常微分方程的数值解法有一个初步了解，同时学会使用MATLAB软件求解微分方程的基本命令，学会建立微分方程方面的数学模型。

这对于学生深入理解微分、积分的数学概念，掌握数学的分析思维方法，熟悉处理大量的工程计算问题的方法是十分必要的。

二、实验过程（实验步骤、记录、数据、分析）1．开启MATLAB软件平台，开启MATLAB编辑窗口；2．根据问题，建立的线性规划模型，并编写求解规划模型的M文件；3．保存文件并运行；4．观察运行结果(数值或图形)，并不断地改变参数设置观察运行结果；5．根据观察到的结果和体会，写出实验报告。

三、实验要求与任务根据实验内容和步骤，完成以下实验，要求写出实验报告（实验目的→问题→数学模型→算法与编程→计算结果→分析、检验和结论）1．求微分方程的解析解，并画出它们的图形。

y '= y + 2 x, y (0) = 1, 0< x <1；程序如下：由y=dsolve('Dy=y+2*x','y(0)=1','x')得出解析解y =-2*x-2+3*exp(x) 建立函数m 文件：function y=myfun4(x)y=-2*x-2+3*exp(x)画图函数为fplot('myfun4',[0,1])图形如下：2．求微分方程⎪⎩⎪⎨⎧====-+]100[0)0(;0)0(01.03t uu u u u 的数值解，要求编写求解程序。

首先建立函数M 文件：function dy=myfun5(t,y)dy=zeros(2,1);dy(1)=y(2);dy(2)=0.1*(y(1).^3)-y(1);输入命令：[T,Y]=ode15s('myfun5',[0,10],[0,1]);plot(T,Y(:,1),'-',T,Y(:,2),'*');图形如下：3．Rossler 微分方程组：⎪⎩⎪⎨⎧-+=+=--=)('''c x z b z ayx y z y x 当固定参数b =2，c =4时，试讨论随参数a 由小到大变化（如 a ∈(0，0.65))而方程解的变化情况，并且画出空间曲线图形，观察空间曲线是否形成混沌状？ 首先建立函数M 文件rossler.m ，在其中用x(1)表示x ，用x(2)表示y ，用x(3)表示z.function r=rossler(t,x)global a;global b;global c;r=[-x(2)-x(3);x(1)+a*x(2);b+x(3)*(x(1)-c)];主程序如下：global a;global b;global c;b=2;c=4;t0=[0,200];for a=0:0.02:0.65[t,x]=ode45('rossler',t0,[0,0,0]);asubplot(1,2,1);plot(t,x(:,1),'r',t,x(:,2),'g',t,x(:,3),'b');title('x(红色),y(绿色),z(篮色)随t变化情况');xlabel('t'); subplot(1,2,2);plot3(x(:,1),x(:,2),x(:,3))title('相图');xlabel('x');ylabel('y');zlabel('z');pauseend当a=0 时，(x,y,z)收敛于（0,0.5,0.5）当a=0.05 时，(x,y,z)仍然收敛，但收敛速度较小。

实验一 常微分方程1. 分别用Euler 法和ode45解下列常微分方程并与解析解比较： (1) ,(0)1,13y x y y x '=+=<<Euler 法：function [t,y]=euler(Fun,tspan,y0,h) t=tspan(1):h:tspan(2); y(1)=y0;for i=1:length(t)-1y(i+1)=y(i)+h.*feval(Fun,t(i),y(i)); end t=t'; y=y';function f=Fun(x,y) % 常微分方程的右端函数 f=x+y;>> [x,y]=euler('Fun',[0,3],1,0.1)>> [x,y] ans =0 1.0000 0.1000 1.1000 0.2000 1.2200 0.3000 1.3620 0.4000 1.5282 0.5000 1.7210 0.6000 1.9431 0.7000 2.1974 0.8000 2.4872 0.9000 2.8159 1.0000 3.1875 1.1000 3.6062 1.2000 4.0769 1.3000 4.6045 1.4000 5.1950 1.5000 5.8545 1.6000 6.5899 1.7000 7.4089 1.8000 8.3198 1.9000 9.3318 2.0000 10.4550 2.1000 11.7005 2.2000 13.0805 2.3000 14.6086 2.4000 16.2995 2.5000 18.1694 2.6000 20.2364 2.7000 22.5200 2.8000 25.0420 2.9000 27.8262 3.0000 30.8988ode45：>> [x,y]=ode45('Fun',[0,3],1) ans =0 1.0000 0.0502 1.0528 0.1005 1.1109 0.1507 1.17460.2010 1.2442 0.2760 1.3596 0.3510 1.4899 0.4260 1.63610.5010 1.7996 0.5760 1.9817 0.6510 2.1838 0.7260 2.4074实验一 常微分方程0.8010 2.6544 0.8760 2.9264 0.9510 3.2254 1.0260 3.55351.1010 3.9131 1.1760 4.3065 1.2510 4.7364 1.3260 5.20561.4010 5.7172 1.4760 6.2744 1.5510 6.8810 1.6260 7.54061.7010 8.2574 1.7760 9.0359 1.8510 9.8808 1.9260 10.79742.0010 11.7912 2.0760 12.8683 2.1510 14.0351 2.2260 15.29862.3010 16.6664 2.3760 18.1466 2.4510 19.7478 2.5260 21.47962.6010 23.3522 2.6760 25.3764 2.7510 27.5641 2.8260 29.92812.9010 32.4820 2.9257 33.3694 2.9505 34.2796 2.9752 35.21343.0000 36.1711解析解：>> y=dsolve('Dy=x+y','y(0)=1','x') y =2*exp(x) - x - 1(2) 20.01()2sin(),(0)0,(0)1,05y y y t y y t ''''-+===<< Euler 法：实验一常微分方程function f=Fun(t,y)% 常微分方程的右端函数f=[y(2);0.01*y(2)^2-2*y(1)+sin(t)];>> [t,y]=euler('Fun',[0,5],[0,1],0.2)ode45：>> [t,y]=ode45('Fun',[0,5],[0,1])t =0 0.0001 0.0001 0.0002 0.0002 0.0005 0.0007 0.0010 0.0012 0.00250.0037 0.0050 0.0062 0.0125 0.0188 0.0251 0.0313 0.0627 0.0941 0.12550.1569 0.2819 0.4069 0.5319 0.6569 0.7819 0.9069 1.0319 1.1569 1.28191.4069 1.5319 1.6569 1.7819 1.90692.0319 2.1569 2.2819 2.4069 2.53192.6569 2.7819 2.90693.0319 3.1569 3.2819 3.4069 3.5319 3.6569 3.78193.90694.0319 4.1569 4.2819 4.4069 4.5319 4.6569 4.7427 4.8285 4.91425.0000y =0 1.0000 0.0001 1.0000 0.0001 1.0000 0.0002 1.0000 0.0002 1.00000.0005 1.0000 0.0007 1.0000 0.0010 1.0000 0.0012 1.0000 0.0025 1.00000.0037 1.0000 0.0050 1.0000 0.0062 1.0000 0.0125 1.0000 0.0188 1.00000.0251 0.9999 0.0313 0.9998 0.0627 0.9987 0.0941 0.9965 0.1253 0.99340.1564 0.9893 0.2786 0.9632 0.3966 0.9220 0.5085 0.8662 0.6126 0.79670.7072 0.7146 0.7908 0.6210 0.8620 0.5176 0.9198 0.4058 0.9632 0.28760.9915 0.1647 1.0043 0.0392 1.0013 -0.0869 0.9826 -0.2117 0.9485 -0.33310.8996 -0.4490 0.8365 -0.5578 0.7605 -0.6577 0.6725 -0.7471 0.5742 -0.8246实验一 常微分方程0.4669 -0.8889 0.3525 -0.9393 0.2327 -0.9748 0.1095 -0.9950 -0.0154 -0.9996-0.1398 -0.9887 -0.2619 -0.9624 -0.3798 -0.9212 -0.4916 -0.8657 -0.5957 -0.7970-0.6904 -0.7161 -0.7742 -0.6242 -0.8460 -0.5228 -0.9046 -0.4134 -0.9491 -0.2978-0.9789 -0.1777 -0.9934 -0.0549 -0.9945 0.0300 -0.9883 0.1146 -0.9748 0.1985-0.9543 0.28092. 求一通过原点的曲线，它在(,)x y 处的切线斜率等于22,0 1.57.x y x +<<若x 上限增为1.58，1.60会发生什么？function f=Fun(x,y) % 常微分方程的右端函数 f=2*x+y.^2;>> [x,y]=ode45('Fun',[0,1.57],0) x =0 0.0393 0.0785 0.1178 0.1570 0.1963 0.2355 0.2748 0.3140 0.3533 0.3925 0.4318 0.4710 0.5103 0.5495 0.5888 0.6280 0.6673 0.7065 0.7458 0.7850 0.8243 0.8635 0.9028 0.9420 0.9813 1.0205 1.0598 1.0990 1.1383 1.1775 1.2168 1.2560 1.2953 1.3345 1.3738 1.4130 1.4248 1.4367 1.4485 1.4604 1.4722 1.4840 1.4959 1.5077 1.5140 1.5203 1.5265 1.5328 1.5376 1.5424 1.5472 1.5519 1.5543 1.5567 1.5591 1.5614 1.5631 1.5647 1.5664 1.5681 1.5685 1.5690 1.5695 1.5700 y =实验一 常微分方程0 0.0015 0.0062 0.0139 0.0247 0.0386 0.0556 0.0758 0.09920.1259 0.1559 0.1895 0.2266 0.2675 0.3124 0.3615 0.4152 0.4738 0.5378 0.6076 0.6841 0.7679 0.8601 0.9620 1.0751 1.2014 1.3434 1.5045 1.6892 1.9037 2.1557 2.4577 2.8282 3.3003 3.9056 4.7317 5.9549 6.4431 7.0116 7.6832 8.4902 9.4821 10.7170 12.3090 14.4551 15.9220 17.7080 19.9390 22.8164 25.6450 29.2282 33.9673 40.5910 44.9434 50.3088 57.1229 66.1087 74.3108 84.7123 98.4901 117.7875 124.9206 132.9699 142.1268 152.641500.20.40.60.81 1.2 1.4 1.6若x 上限增为1.58，1.60,则超出运算的范围，发生溢出。