数学实验 常微分方程
常微分方程数值解实验报告
常微分方程数值解实验报告学院：数学与信息科学专业：信息与计算科学姓名：郑思义学号：201216524课程：常微分方程数值解实验一：常微分方程的数值解法1、分别用Euler 法、改进的Euler 法（预报校正格式）和S —K 法求解初值问题。
（h=0.1）并与真解作比较。
⎩⎨⎧=++-=10(1y'）y x y 1.1实验代码：%欧拉法function [x,y]=naeuler(dyfun,xspan,y0,h)%dyfun 是常微分方程，xspan 是x 的取值范围，y0是初值，h 是步长 x=xspan(1):h:xspan(2); y(1)=y0;for n=1:length(x)-1y(n+1)=y(n)+h*feval(dyfun,x(n),y(n)); end%改进的欧拉法function [x,m,y]=naeuler2(dyfun,xspan,y0,h)%dyfun 是常微分方程，xspan 是x 的取值范围，y0是初值，h 是步长。
%返回值x 为x 取值，m 为预报解，y 为校正解 x=xspan(1):h:xspan(2); y(1)=y0;m=zeros(length(x)-1,1); for n=1:length(x)-1 k1=feval(dyfun,x(n),y(n)); y(n+1)=y(n)+h*k1; m(n)=y(n+1);k2=feval(dyfun,x(n+1),y(n+1)); y(n+1)=y(n)+h*(k1+k2)/2; end%四阶S —K 法function [x,y]=rk(dyfun,xspan,y0,h)%dyfun 是常微分方程，xspan 是x 的取值范围，y0是初值，h 是步长。
x=xspan(1):h:xspan(2); y(1)=y0;for n=1:length(x)-1 k1=feval(dyfun,x(n),y(n));k2=feval(dyfun,x(n)+h/2,y(n)+(h*k1)/2); k3=feval(dyfun,x(n)+h/2,y(n)+(h*k2)/2); k4=feval(dyfun,x(n)+h,y(n)+h*k3);y(n+1)=y(n)+(h/6)*(k1+2*k2+2*k3+k4);end%主程序x=[0:0.1:1];y=exp(-x)+x;dyfun=inline('-y+x+1');[x1,y1]=naeuler(dyfun,[0,1],1,0.1);[x2,m,y2]=naeuler2(dyfun,[0,1],1,0.1);[x3,y3]=rk(dyfun,[0,1],1,0.1);plot(x,y,'r',x1,y1,'+',x2,y2,'*',x3,y3,'o');xlabel('x');ylabel('y');legend('y为真解','y1为欧拉解','y2为改进欧拉解','y3为S—K解','Location','NorthWest');1.2实验结果：x 真解y 欧拉解y1 预报值m 校正值y2 S—K解y30.0 1.0000 1.0000 1.0000 1.00000.1 1.0048 1.0000 1.0000 1.0050 1.00480.2 1.0187 1.0100 1.0145 1.0190 1.01870.3 1.0408 1.0290 1.0371 1.0412 1.04080.4 1.0703 1.0561 1.0671 1.0708 1.07030.5 1.1065 1.0905 1.1037 1.1071 1.10650.6 1.1488 1.1314 1.1464 1.1494 1.14880.7 1.1966 1.1783 1.1945 1.1972 1.19660.8 1.2493 1.2305 1.2475 1.2500 1.24930.9 1.3066 1.2874 1.3050 1.3072 1.30661.0 1.3679 1.3487 1.3665 1.3685 1.36792、选取一种理论上收敛但是不稳定的算法对问题1进行计算，并与真解作比较。
实验报告——常微分方程的数值解法
实验报告实验项目名称常微分方程的数值解法实验室数学实验室所属课程名称微分方程数值解实验类型上机实验实验日期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\）的变化。
常微分方程数值解实验报告
常微分方程数值解实验报告实验报告：常微分方程数值解1.引言常微分方程（Ordinary Differential Equations, ODEs）是数学领域中一个重要的研究对象，涉及到许多自然科学和工程技术领域的问题。
解常微分方程的数值方法是一种求解差分方程的方法，通过计算机找到方程的近似解，对于模拟和预测连续过程非常有用。
本实验旨在通过数值解法，验证和应用常微分方程的解，并比较不同数值方法的精度和效率。
2.实验目的2.1理解常微分方程的基本概念和数值解法；2.2掌握将常微分方程转化为数值求解问题的基本方法；2.3运用数值解法求解常微分方程；2.4比较不同数值解法的精度和效率。
3.实验内容3.1 欧拉方法（Euler Method）给定一个一阶常微分方程dy/dx=f(x,y)，通过将其离散为差分形式，欧拉方法可以通过以下递推公式来求解：y_{n+1}=y_n+h*f(x_n,y_n)其中，h为步长，x_n和y_n为当前的x和y值。
3.2 改进的欧拉方法（Improved Euler Method）改进的欧拉方法使用欧拉方法的斜率的平均值来估计每一步中的斜率。
具体公式如下：k1=f(x_n,y_n)k2=f(x_n+h,y_n+h*k1)y_{n+1}=y_n+h*((k1+k2)/2)3.3 二阶龙格-库塔法（Second-order Runge-Kutta Method）二阶龙格-库塔法通过计算每个步骤中的两个斜率来估计每个步长中的斜率。
具体公式如下：k1=f(x_n,y_n)k2=f(x_n+h/2,y_n+(h/2)*k1)y_{n+1}=y_n+h*k24.实验步骤4.1选取常微分方程，并将其转化为数值求解问题的形式；4.2根据给定的初始条件和步长，使用欧拉方法、改进的欧拉方法和二阶龙格-库塔法求解该方程；4.3比较三种方法的数值解与理论解的差异，并分析其精度和效率；4.4尝试不同的步长，观察相应的数值解的变化。
常微分方程数值解实验
多步法，Gear’s反向
数值积分，精度中等
若ode45失效时，
可尝试使用
ode23s
刚性
一步法，2阶Rosebrock算法，
低精度。
当精度较低时，
计算时间比ode15s短
odefx为显式常微分方程 中的 ，t为求解区间，要获得问题在其他指定点 上的解，则令t=[t0,t1,t2,…]（要求 单调），y0初始条件。
MATLAB 中有几个专门用于求解常微分方程的函数,它们的设计思想基于Runge-Kutta方法,基本设计思想为：从改进的欧拉方法比欧拉方法精度高的缘由着手,如果在区间[ x1, xi+1]多取几个点的斜率值,然后求取平均值,则可以构造出精度更高的计算方法。 这些函数主要包括：ode45、ode23、ode15s、ode113、ode23s、ode23t、ode23tb. 其中最常用的是函数ode45,该函数采用变步长四阶五阶Runge-Kutta法求数值解,并采用自适应变步长的求解方法。ode23采用二阶三阶Runge-Kutta法求数值解,与ode45类似,只是精度低一些。ode15s用来求刚性方程组。
43
4月22日
588
666
28
46
4月23日
693
782
35
55
4月24日
774
863
39
64
4月25日
877
954
42
73
4月26日
988
1093
48
76
4月27日
1114
1255
56
78
4月28日
1199
1275
59
78
4月29日
常微分方程实验报告
常微分方程实验报告《常微分方程》综合性实验实验报告实验班级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．了解常微分⽅程的基本概念。
2．了解常微分⽅程的解析解。
3．了解常微分⽅程的数值解。
4．学习掌握MATLAB 软件有关的命令。
【实验内容】如右图所⽰，⼀根长l 的⽆弹性细线，⼀段固定，另⼀端悬挂⼀个质量为m 的⼩球，在重⼒的作⽤下⼩球处于垂直的平衡位置。
若使⼩球偏离平衡位置⼀个⾓度θ，让它⾃由，它就会沿圆弧摆动。
在不考虑空⽓阻⼒的情况下，⼩球会做⼀定周期的简谐运动。
利⽤⽜顿第⼆定律得到如下的微分⽅程0)0(',)0(,sin "0===θθθθθmg ml问该微分⽅程是线性的还是⾮线性的？是否存在解析解？如果不存在解析解，能否求出其近似解？【实验准备】1．微分⽅程的概念未知的函数以及它的某些阶的导数连同⾃变量都由⼀已知⽅程联系在⼀起的⽅程称为微分⽅程。
如果未知函数是⼀元函数，称为常微分⽅程。
常微分⽅程的⼀般形式为0),,",',,()(=n y y y y t F如果未知函数是多元函数，成为偏微分⽅程。
联系⼀些未知函数的⼀组微分⽅程组称为微分⽅程组。
微分⽅程中出现的未知函数的导数的最⾼阶解数称为微分⽅程的阶。
若⽅程中未知函数及其各阶导数都是⼀次的，称为线性常微分⽅程，⼀般表⽰为)()(')()(1)1(1)(t b y t a y t a y t a y n n n n =++++--若上式中的系数n i t a i ,,2,1),( =均与t ⽆关，称之为常系数或定常、⾃治、时不变的。
2．常微分⽅程的解析解有些微分⽅程可直接通过积分求解.例如,⼀解常系数常微分⽅程1+=y dtdy可化为dt y dy=+1,两边积分可得通解为1-=t ce y .其中c 为任意常数.有些常微分⽅程可⽤⼀些技巧,如分离变量法,积分因⼦法,常数变异法,降阶法等可化为可积分的⽅程⽽求得解析解(显式解).线性常微分⽅程的解满⾜叠加原理,从⽽他们的求解可归结为求⼀个特解和相应齐次微分⽅程的通解.⼀阶变系数线性微分⽅程总可⽤这⼀思路求得显式解。
大学数学实验四_常微分方程数值解
2 小型火箭初始质量为 1400kg，其中包括 1080kg 燃料。火箭竖直向上发射时，燃料燃烧率 为 18kg/s，由此产生 32000N 推力，火箭引擎在燃料用尽时关闭。设火箭上升时空气阻力正 比于速度的平方，比例系数为 0.4kg/m，求引擎关闭瞬间火箭的高度、速度、加速度，及火 箭到达最高点时的高度和加速度，并画出高度、速度、加速度随时间变化的图形。
6
输出结果为
故，当 t =60s 时，火箭的速度、高度、加速度为
v (m/s)
向前欧拉公式
266.94
经典龙格-库塔公式
266.94
ode45(@rocket,t,v)
266.93
h (m) 12131 12131 12131
从上表可以看出，用三种不同方法求得的结果相差不大。
60s 后，火箭作竖直上抛运动，加速度为
4
输出的结果为
输出的火箭的 h – t 图线如下：
经典龙格-库塔公式： n=100000; h=60/n; t=0:h:60; v(1)=0; z(1)=0; g(1)=90/7; for i=1:n
a = v(i); s = z(i); x1 = h*i; k1 = (32000-0.4*a.^2)/(1400-18*x1)-10; x1 = x1 + h/2; b = a + h/2*k1; k2 = (32000-0.4*b.^2)/(1400-18*x1)-10; c = a + h/2*k2; k3 = (32000-0.4*c.^2)/(1400-18*x1)-10; x1 = x1 + h/2;
n = 500 267.0539 n = 100000 266.9432
数学实验基础实验报告常微分方程
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,则超出运算的范围，发生溢出。
数学实验-常微分方程
实验六 常微分方程的Matlab 解法一、实验目的1． 了解常微分方程的解析解。
2． 了解常微分方程的数值解。
3． 学习掌握MATLAB 软件有关的命令。
二、实验内容一根长l 的无弹性细线，一段固定，另一端悬挂一个质量为m 的小球，在重力的作用下小球处于垂直的平衡位置。
若使小球偏离平衡位置一个角度θ，让它自由，它就会沿圆弧摆动。
在不考虑空气阻力的情况下，小球会做一定周期的简谐运动。
利用牛顿第二定律得到如 下的微分方程0)0(',)0(,sin "0===θθθθθmg ml问该微分方程是线性的还是非线性的？是否存在解析解？如果不存在解析解，能否求出其近似解？三、实验准备MATLAB 中主要用dsolve 求符号解析解，ode45,ode23,ode15s 求数值解。
ode45是最常用的求解微分方程数值解的命令，对于刚性方程组不宜采用。
ode23与ode45类似，只是精度低一些。
ode12s 用来求解刚性方程组，是用格式同ode45。
可以用help dsolve, help ode45查阅有关这些命令的详细信息.四、实验方法与步骤练习1 求下列微分方程的解析解 （1）b ay y +='（2）1)0(',0)0(,)2sin(''==-=y y y x y （3）1)0(',1)0(',','==-=+=g f f g g g f f 方程（1）求解的MATLAB 代码为：clear;s=dsolve('Dy=a*y+b')结果为s =-b/a+exp(a*t)*C1方程（2）求解的MATLAB 代码为：clear;s=dsolve('D2y=sin(2*x)-y','y(0)=0','Dy(0)=1','x') simplify(s) %以最简形式显示s结果为s =(-1/6*cos(3*x)-1/2*cos(x))*sin(x)+(-1/2*sin(x)+1/6*sin(3*x))*cos(x)+5/3*sin(x) ans =-2/3*sin(x)*cos(x)+5/3*sin(x) 方程（3）求解的MATLAB 代码为：clear;s=dsolve('Df=f+g','Dg=g-f','f(0)=1','g(0)=1') simplify(s.f) %s 是一个结构 simplify(s.g)结果为ans =exp(t)*cos(t)+exp(t)*sin(t) ans =-exp(t)*sin(t)+exp(t)*cos(t) 练习２ 求解微分方程,1)0(,1'=++-=y t y y先求解析解，再求数值解，并进行比较。
数学实验基础 实验报告(1)常微分方程
实验一 常微分方程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,则超出运算的范围，发生溢出。
常微分方程实验
常微分方程实验
1.1 实验目的
理解常微分方程及常微分方程组的概念，掌握 MATLAB 命 令求一阶、二阶常微分方程及常微分方程组的通解与满足一定初 始或边界条件的特解。
1.2 常微分方程和常微分方程组的求解
MATLAB 中，若 y 是因变量，用“Dny”表示 y 的 n 阶导， 如 Dy 表示 y′ ，D2y表示 y″ ，Dy (0) = 1表示 y′ (0) = 1。常微分方 程的求解可以调用 dsolve 函数来实现，其一般格式为
dsolve ( ′ eq ′ , ′ c ′ , ′ x ′ )：求微分方程 eq 在初始条件 c 下的 特解，参数 x 为微分方程中的自变量，省略参数 c 时，则求方程 的通解。
dsolve ( ′ eq1 ′ , ′ eq2 ′ , … ,′ eqn ′ ,′ c1 ′ , ′ c2 ′ , … , ′ cn ′ , ′ x ′ )： 求微分方程组 eq1，eq2，…，eqn 在初始条件c1，c2，…，cn下 的特解，参数 x 为微分方程组中的自变量，省略初始条件，则求 方程组的通解。
%求(2)式常微分方程的特解
输出：ans = 2 / ( exp ( -x ) -2 * exp ( x ))
dsolve ( ′ D2y -6 * Dy + 13 * y = 0 ′ , ′ x ′ )
%求(3)式二
阶常微分方程的通解
输出：ans = C1 * exp ( 3 * x ) * sin ( 2 * x ) + C2 * exp ( 3 * x )
高等数学
解 在 MATLAB 指令窗口输入命令如下：
syms x y t；
dsolve ( ′ Dy + y / x = sin ( x ) / x ′ , ′ x ′ )
数学实验第二次作业任务常微分方程数值求解
数学实验第⼆次作业任务常微分⽅程数值求解实验4常微分⽅程数值解实验⽬的：1.练习数值积分的计算；2.掌握⽤MATLAB软件求微分⽅程初值问题数值解的⽅法；3.通过实例学习⽤微分⽅程模型解决简化的实际问题；4.了解欧拉⽅法和龙格——库塔⽅法的基本思想和计算公式，及稳定性等概念。
实验内容：3.⼩型⽕箭初始质量为1400kg，其中包括1080kg燃料，⽕箭竖直向上发射是燃料燃烧率为18kg/s，由此产⽣32000N的推⼒，⽕箭引擎在燃料⽤尽时关闭。
设⽕箭上升是空⽓阻⼒正⽐于速度的平⽅，⽐例系数为0.4kg/m，求引擎关闭瞬间⽕箭的⾼度，速度，加速度，及⽕箭到达最⾼点是的⾼度，速度和加速度，并画出⾼度，速度，加速度随时间变化的图形。
解答如下：这是⼀个典型的⽜顿第⼆定律问题，分析⽕箭受⼒情况；先规定向上受⼒为正数建⽴数学模型：A燃料未燃尽前，在任意时刻（t<60s）⽕箭受到向上的-F=32000N，向下的重⼒G=mg，g=9.8，向下的阻⼒f=kv^2, k=0.4, v表⽰此时⽕箭速度；此时⽕箭收到的合⼒为F1=（F-mg-f）；⽕箭的初始质量为1400kg，燃料燃烧率为-18kg/s;此刻⽕箭质量为m=1400-18*t根据⽜顿第⼆定律知，加速度a=F1/m=(F-mg-f)/(m-r*t)=(32000-(0.4.*v.^2)-9.8.*(1400-18.*t))由此可利⽤龙格-库塔⽅法来实现，程序实现如下Function [dx]=rocket[t,x] %建⽴名为rocket的⽅程m=1400;k=0.4;r=-18;g=9.8; %给出题⽬提供的常数值dx=[x(2);(32000-(k*x(2)^2)-g*(m+r*t))/(m+r*t)];%以向量的形式建⽴⽅程[a]=(32000-(k*x(2)^2)-g*(m+r*t))/(m+r*t); %给出a的表达式End;ts=0:60; %根据题⽬给定燃烧率计算出燃料燃尽的时间，确定终点x0=[0,0]; %输⼊x的初始值[t,x]=ode15s(@rocket,ts,x0); %调⽤ode15s计算[t,x];h=x(:,1);v=x(:,2);plot(t,x(:,1)),grid; %绘出⽕箭⾼度与时间的关系曲线title('h-t');xlabel('t/s');ylabel('h/m'),pause;plot(t,x(:,2)),grid ; %绘出⽕箭速度与时间的曲线关系title('v-t');xlabel('t/s');ylabel('v/m/s'),pause;a=(32000-(0.4.*v.^2)-9.8.*(1400-18.*t))/(1400-18.*t);plot(t,a),grid; %绘出⽕箭加速度与时间的曲线关系title('a-t'); xlabel('t/s'),ylabel('a/m^2/s'),pause⽕箭⾼度随时间变化的曲线⽕箭速度随时间变化的曲线⽕箭加速度随时间变化的曲线数据过多，故截取部分如下第⼀列为时间，第⼆列为⽕箭⾼度，第三列为⽕箭速度由此可以，在t=60s时，即⽕箭燃料燃尽瞬间，引擎关闭瞬间，⽕箭将到达12912m的⾼度，速度为267,29m，加速度a=0.9m/s^2B燃料燃尽之后，与A 类似，分析受⼒如下⽕箭受到向上的F=0向下的重⼒G=mg，g=9.8，向下的阻⼒f=kv^2, k=0.4, v表⽰此时⽕箭速度；此时⽕箭收到的合⼒为F2=（-mg-f）；⽕箭的初始质量为320kg，恒定根据⽜顿第⼆定律,加速度a=F2/m=-g-0.4v^2/320;程序实现如下function [ dx ] = rocket2( t,x ) %建⽴以rocket2为名的函数dx=[x(2);-9.8-0.4.*x(2).^2/320]; %以向量的形式建⽴⽅程ts=60:120; %给出初始时刻，估计终点时刻x0=[12190,267.26]; %给出x初始值[t,x]=ode15s(@rocket2,ts,x0); %调⽤ode15s计算[t,x]plot(t,x(:,1)),grid; %绘出⽕箭⾼度随时间变化的曲线title('h-t');xlabel('t/s'),ylabel('h/m'),pause;plot(t,x(:,2)),grid; %绘出⽕箭速度随时间的变化曲线title('v-t');xlabel('t/s'),ylabel('v/m/s'),pause;v=x(:,2);a=-9.8-0.4*v.^2/320; %给出加速度的具体表达式plot(t,a),grid; %绘出⽕箭加速度随时间变化的曲线title('a-t'); xlabel('t/s'),ylabel('a/m^2/s'),pause得到的曲线图形如下⽕箭⾼度随时间的变化曲线从图中可以⼤致看出，最⾼点在13km左右，⽕箭速度随时间的变化曲线加速度随时间变化曲线如下数据表格⼤致如下从图表中可以看出，在71s左右速度到达0，即此时到达最⾼处，⾼度为13117m加速度为-9.8m/m/s^2;本题总结：这道题是典型的物理⽜顿⼒学的题⽬，通过受⼒的正确分析，可以知道，以[h,v]为向量建⽴微分⽅程即可求解，h的微分是速度v，速度v的微分是加速度a解题过程中存在的难点是：取值步长不太容易确定，⽽且是哪种算法不确定，先⽤ode15s速度较快，ode23s速度差不太多，其他两种速度较慢，等待时间较长5．⼀只⼩船渡过宽为d 的河流，⽬标是起点A 正对着的另⼀岸B 点。
计算方法--常微分方程求解实验
实验五 常微分方程求解实验一、 实验目的通过本实验学会对给定初值我呢他，用欧拉法、改进欧拉法、四阶龙格-库塔法求数值解和误差，并比较优缺点.对给定刚性微分方程，求其数值解，并与精确解比较，分析计算结果.二、 实验题目1. 解初值问题各种方法比较 实验题目：给定初值问题⎪⎩⎪⎨⎧=≤<+=,0)1(,21,e d d y x x xyx y x 精确解为)e -e (xx y =，按 （1） 欧拉法，步长;1.0,025.0==h h （2） 改进欧拉法，步长;01.0,05.0==h h （3） 四阶标准龙格-库塔法，步长;1.0=h求在节点)10,,2,1(1.01 =+=k k x k 处的数值解及误差，比较各方法的优缺点. 2. 刚性方程计算实验题目：给定刚性微分方程⎪⎩⎪⎨⎧=≤<-+-=-,2)0(,50,600e 8.1199600d d 1.0y x y xyx 其精确解为12e e )(-0.1x -600x-+=x y .任选取一显示方法，取不同的步长求解，并分析计算结果.三、 实验原理1.欧拉格式由数值微分的向前差商公式可以解决初值问题（6.1）()()⎩⎨⎧=≤≤=00',,,y x y b x a y x f y 中的导数的数值计算问题：()()()()(),'1h x y x y h x y h x y x y n n n n -=-+≈+由此可得()()().'1n n n x hy x y x y +≈+（6.1）实际上给出()()()()()().,','n n n x y x f x y x y x f x y =⇒=于是有()()()().,1n n n n x y x hf x y x y +≈+再由()()11,++≈≈n n n n x y y x y y 得().,1,0,,111 =+=+++n y x hf y y n n n n （6.2） 递推公式（6.2）称为欧拉格式。
实验4：常微分方程实验
撰写人姓名： 撰写时间： 审查人姓名：实 验 全 过 程 记 录实验名称常微分方程实验时间地点813姓 名 张强 学 号 1012010817同实验者学 号一、实验目的1、熟练掌握基本的常微分方程的解析求解方法；2、熟练掌握基本的常微分方程的数值求解方法；3、能利用MATLAB 较熟练地求解简单的常微分方程应用问题，并根据需要进行相应 讨论。
二、实验内容：1、练习利用MATLAB 解析求解基本的常微分方程问题；2、练习利用MATLAB 数值求解基本的常微分方程问题；3、练习利用MATLAB 求解简单的常微分方程应用问题。
三、实验用仪器设备及材料软件需求：操作系统：Windows XP 或更新的版本； 实用数学软件：MATLAB 7.0或更新的版本。
硬件需求：Pentium IV 450以上的CPU 处理器、512MB 以上的内存、5000MB 的自由硬盘空间、 CD-ROM 驱动器、打印机、打印纸等。
四、实验原理:高等数学、数值分析等相关理论五、实验步骤：1、 一质量为m 的质点作直线运动，从速度等于0的时刻起，有一个和时间成正比的力（比例系数为k1）作用在上面，此外又受到介质的阻力（阻力与速度成正比，比例系数k2），求质点的速度与时间的关系。
>> syms m k1 k2 v t;>> dsolve('k1*t-k2*v=m*Dv') ans =exp(-k2/m*t)*C1+k1*(k2*t-m)/k2^22、给定如下常微分方程，利用dsolve 函数，求其解析解。
⑴ 25743x x x x t t +-+=--；>> syms t>> x=dsolve('5*D3x+7*D2x-4*Dx+x=-t^2-3*t') x =-30-11*t-t^2+C1*exp(1/1425720*(-47524*(10484+60*7509^(1/2))^(1/3)-2621*(10484+60*7509^(1/2))^(2/3)+15*(10484+60*7509^(1/2))^(2/3)*7509^(1/2)-665336)*t )+C2*exp(-1/2851440*(-47524*(10484+60*3^(1/2)*2503^(1/2))^(1/3)-2621*(10484+60*3^(1/2)*2503^(1/2))^(2/3)+1330672+15*(10484+60*3^(1/2)*2503^(1/2))^(2/3)*3^(1/2)*2503^(1/2))*t)*cos(1/2851440*(10484+60*3^(1/2)*2503^(1/2))^(1/3)*(2621*3^(1/2)*(10484+60*3^(1/2)*2503^(1/2))^(1/3)-45*(10484+60*3^(1/2)*2503^(1/2))^(1/3)*2503^(1/2)-47524*3^(1/2))*t)+C3*exp(-1/2851440*(-47524*(10484+60*3^(1/2)*2503^(1/2))^(1/3)-2621*(10484+60*3^(1/2)*2503^(1/2))^(2/3)+1330672+15*(10484+60*3^(1/2)*2503^(1/2))^(2/3)*3^(1/2)*2503^(1/2))*t)*sin(1/2851440*(10484+60*3^(1/2)*2503^(1/2))^(1/3)*(2621*3^(1/2)*(10484+60*3^(1/2)*2503^(1/2))^(1/3)-45*(10484+60*3^(1/2)*2503^(1/2))^(1/3)*2503^(1/2)-47524*3^(1/2))*t)⑵ ⎪⎩⎪⎨⎧='=-=''==012cos 00x x y y y x y 。
数学建模-- 常微分方程数值解及实验
注意:
1、在解n个未知函数的方程组时，x0和x均为n维向量， m-文件中的待解方程组应以x的分量形式写成. 2、使用Matlab软件求数值解时，高阶微分方程必须 等价地变换成一阶微分方程组.
例
d 2x dx 2 x 0 2 ( x 1) dt dt x (0 ) 2; x '(0 ) 0
实际应用时，与欧拉公式结合使用：
0 y i( 1) y i hf ( x i , y i ) h ( k 1) (k ) y i 1 y i [ f ( x i , y i ) f ( x i 1 , y i 1 )] k 0 ,1, 2 , 2
•欧拉法是一阶公式，改进的欧拉法是二阶公式。
•龙格-库塔法有二阶公式和四阶公式。 •线性多步法有四阶阿达姆斯外插公方程的数值解
[t，x]=solver（’f’,ts,x0,options）
自变 量值 函数 值
ode45 ode23 ode113 ode15s ode23s
1、建立m-文件vdp1.m如下： function dy=vdp (t,y) dy=zeros(2,1); dy(1)=y(2); dy(2)=(1-y(1)^2)*y(2)-y(1); 2、取t0=0，tf=20，输入命令： [T,Y]=ode15s('vdp1',[0 20],[2 0]); plot(T,Y(:,1),'-')
y ( x i 1 ) y ( x i )
xi 1
f ( t , y ( t )) dt
x i 1 x i 2
xi
[ f ( x i , y ( x i )) f ( x i 1 , y ( x i 1 ))]
常微分方程的数值解法实验报告
常微分方程的数值解法实验报告实验报告：常微分方程的数值解法摘要：常微分方程（ODE）是描述动力学系统中物理量随时间变化的数学方程，广泛应用于自然科学和工程领域。
然而，对于一些复杂的非线性ODE，很难找到解析解。
因此，我们需要数值解法来求解这些方程。
本实验报告将介绍四种常见的常微分方程数值解法，分别是欧拉法、改进的欧拉法、四阶龙格-库塔法和自适应步长的龙格-库塔法，并通过数值实验比较它们的精度和效率。
1.引言在实际问题中，许多物理量的变化规律可以由常微分方程描述。
然而，对于复杂的非线性ODE，很难找到解析解。
因此，为了解决这类问题，我们需要借助数值方法来求解。
2.方法本实验采用四种常见的常微分方程数值解法：欧拉法、改进的欧拉法、四阶龙格-库塔法和自适应步长的龙格-库塔法。
（1）欧拉法是最简单的数值解法，通过将微分方程转化为差分方程，使用离散的步长来近似微分方程。
（2）改进的欧拉法在欧拉法的基础上进行了改进，使用预估-校正的方法来提高精度。
（3）四阶龙格-库塔法是一种经典的数值解法，通过利用不同步长处的斜率来近似微分方程，具有较高的精度。
（4）自适应步长的龙格-库塔法是在四阶龙格-库塔法的基础上改进而来的，根据步长的大小自适应地选择不同的步长，同时保证精度和效率。
3.实验设计为了比较这四种数值解法的精度和效率，我们设计了两个实验。
实验一是求解一阶常微分方程：dy/dx = -2x，初始条件y(0) = 1，解析解为y = 1 - x^2、实验二是求解二阶常微分方程：d^2y/dx^2 + y = 0，初始条件y(0) = 0，dy/dx(0) = 1，解析解为y = sin(x)。
4.结果与分析实验一中，比较四种数值解法在不同步长下的近似解和解析解，计算其误差。
实验结果表明，四阶龙格-库塔法和自适应步长的龙格-库塔法具有较高的精度，而欧拉法和改进的欧拉法的精度较低。
实验二中，我们比较四种数值解法在不同步长下的近似解和解析解，并计算其误差。
常微分方程的求解与定性分析实验报告
常微分方程的求解与定性分析实验报告 一、实验综述1、实验目的及要求● 归纳和学习求解常微分方程(组)的基本原理和方法；● 掌握解析、数值解法，并学会用图形观察解的形态和进行解的定性分析； ● 熟悉MATLAB 软件关于微分方程求解的各种命令；● 通过范例学习建立微分方程方面的数学模型以及求解全过程；● 通过该实验的学习，使学生掌握微分方程(组)求解方法（解析法、欧拉法、梯度法、改进欧拉法等），对常微分方程的数值解法有一个初步了解，同时学会使用MATLAB 软件求解微分方程的基本命令，学会建立微分方程方面的数学模型。
这对于学生深入理解微分、积分的数学概念，掌握数学的分析思维方法，熟悉处理大量的工程计算问题的方法是十分必要的。
2、实验仪器、设备或软件电脑 、matlab7.0二、实验过程（实验步骤、记录、数据、分析）实验内容：根据实验内容和步骤，完成以下实验，要求写出实验报告（实验目的→问题→数学模型→算法与编程→计算结果→分析、检验和结论）1．求微分方程的解析解，并画出它们的图形。
y '= y + 2 x, y (0) = 1, 0< x <1；m=dsolve('Dy=y+2*x','y(0)=1','x')ezplot(m,[0 1])m =3*exp(x) - 2*x – 21．求微分方程⎪⎩⎪⎨⎧====-+]100[0)0(;0)0(01.03t uu u u u &&& 的数值解，要求编写求解程序。
function dy=vdp1000(t,y)dy=zeros(2,1);dy(1)=y(2);dy(2)=-y(1)+0.1*y(1)^3;[T,Y]=ode15s('vdp1000',[0 10],[0 0]);plot(T,Y(:,1),'-')3．Rossler 微分方程组：当固定参数b =2，c =4时，试讨论随参数a 由小到大变化（如 a ∈(0，0.65))而方程解的变化情况，并且画出空间曲线图形，观察空间曲线是否形成混沌状？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.1:0.6[t,x]=ode45('rossler',t0,[0,0,0]);subplot(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:a=0.1:a=0.2:a=0.3:a=0.4:a=0.5:结果分析：从图像可以看出，当a=0时，微分方程的解（x,y,z）收敛与（0，0.5，0.5）；当a=0.1时，（x,y,z）仍收敛与（0，0.5，0.5），只是收敛速度减慢；当a=0.2时，（x,y,z）已发散，周期性变化；随着a的增大，(x,y,z)接近其极限环的速度加快，空间曲线成混沌状。
常微分方程实验
常微分方程实验实验4 研究如下的Lorenz 方程⎪⎩⎪⎨⎧-=--=-=bz xy z xz y rx y y x a x '')('如果你是第一次接触这个方程，相信你一定觉得它并无特别之处。
千万不要小看了这个简单的问题！Lorenz 方程具有重要的实际背景，它所揭示的现象---著名的Lorenz 吸引子，具有重大的科学价值。
本实验就将引导你进入混沌的殿堂。
吸引子实验的背景前面研究的离散动力系统，可以看成是在持续的离散时间间隔上，一个数 量值的变化关系。
如果时间的间隔可以趋于零，那么这种变化关系通常就可以 用微分方程描述。
这样，很自然地就可以把一个(时间独立的)微分方程作为一 个连续的动力系统进行研究。
这里我们最关注的是当时间趋于无穷时，对于某 个集合上的初始值，微分方程解的性态。
设 D 是 n R 中的一个区域，再设 f :Dn R 是一个光滑函数。
考虑微分方程'()()x t f x = (*) 吸引子是描述连续动力系统长期性质的重要概念。
区域 D 的一个闭子集 F 叫 做包含 F 的吸引域 V 的吸引子，如果对开集 V 中的所有初始点 x(0), 通过 x(0)的轨道 x (t)随着 t 的无限增大而接近 F 。
当然，应要求 F 是不变的， 即 如果 x(0)是 F 中的点，那么对0t < t <+∞, ()x t 都在 F 中，这就意味着 F 是一些轨道的并集。
同时我们还要求 在下面的意义下是最小的，即存在一点 (0)x , 使得 ()x t 在 F 中稠密。
具有复杂吸引子连续动力系统最简单又引人入胜的例子，就是Lorenz 方程。
该系统表面上并无任何神秘之处，但它有很强的实际背景，它是长期气象预报的 一个简化模型。
⎪⎩⎪⎨⎧-=--=-=bz xy z xz y rx y y x a x '')('在 (x,y,z) 空间中考虑我们的问题。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
实验六 常微分方程的Matlab 解法一、实验目的1． 了解常微分方程的解析解。
2． 了解常微分方程的数值解。
3． 学习掌握MATLAB 软件有关的命令。
二、实验内容一根长l 的无弹性细线，一段固定，另一端悬挂一个质量为m 的小球，在重力的作用下小球处于垂直的平衡位置。
若使小球偏离平衡位置一个角度θ，让它自由，它就会沿圆弧摆动。
在不考虑空气阻力的情况下，小球会做一定周期的简谐运动。
利用牛顿第二定律得到如 下的微分方程0)0(',)0(,sin "0===θθθθθmg ml问该微分方程是线性的还是非线性的？是否存在解析解？如果不存在解析解，能否求出其近似解？三、实验准备MATLAB 中主要用dsolve 求符号解析解，ode45,ode23,ode15s 求数值解。
ode45是最常用的求解微分方程数值解的命令，对于刚性方程组不宜采用。
ode23与ode45类似，只是精度低一些。
ode12s 用来求解刚性方程组，是用格式同ode45。
可以用help dsolve, help ode45查阅有关这些命令的详细信息.四、实验方法与步骤练习1 求下列微分方程的解析解 （1）b ay y +='（2）1)0(',0)0(,)2sin(''==-=y y y x y （3）1)0(',1)0(',','==-=+=g f f g g g f f 方程（1）求解的MATLAB 代码为：clear;s=dsolve('Dy=a*y+b')结果为s =-b/a+exp(a*t)*C1方程（2）求解的MATLAB 代码为：clear;s=dsolve('D2y=sin(2*x)-y','y(0)=0','Dy(0)=1','x') simplify(s) %以最简形式显示s结果为s =(-1/6*cos(3*x)-1/2*cos(x))*sin(x)+(-1/2*sin(x)+1/6*sin(3*x))*cos(x)+5/3*sin(x) ans =-2/3*sin(x)*cos(x)+5/3*sin(x) 方程（3）求解的MATLAB 代码为：clear;s=dsolve('Df=f+g','Dg=g-f','f(0)=1','g(0)=1') simplify(s.f) %s 是一个结构 simplify(s.g)结果为ans =exp(t)*cos(t)+exp(t)*sin(t) ans =-exp(t)*sin(t)+exp(t)*cos(t) 练习２ 求解微分方程,1)0(,1'=++-=y t y y先求解析解，再求数值解，并进行比较。
由clear;s=dsolve('Dy=-y+t+1','y(0)=1','t') simplify(s)可得解析解为te t y -+=。
下面再求其数值解，先编写M 文件fun8.m%M 函数fun8.mfunction f=fun8(t,y) f=-y+t+1;再用命令clear; close; t=0:0.1:1;y=t+exp(-t); plot(t,y); %化解析解的图形hold on; %保留已经画好的图形，如果下面再画图，两个图形和并在一起 [t,y]=ode45('fun8',[0,1],1);plot(t,y,'ro'); %画数值解图形，用红色小圈画 xlabel('t'),ylabel('y')结果见图6.7.1图6.7.1 解析解与数值解由图7.1可见，解析解和数值解吻合得很好。
下面我们讨论实验引例中的单摆问题.练习3 求方程0)0(',)0(,sin "0===θθθθθmg ml的数值解.不妨取15)0(,8.9,1===θg l .则上面方程可化为0)0(',15)0(,sin 8.9"===θθθθ先看看有没有解析解.运行MATLAB 代码clear;s=dsolve('D2y=9.8*sin(y)','y(0)=15','Dy(0)=0','t') simplify(s)知原方程没有解析解.下面求数值解.令',21θθ==y y 可将原方程化为如下方程组⎪⎩⎪⎨⎧====0)0(,15)0()sin(8.9''211221y y y y y y建立M 文件fun9.m 如下%M 文件fun9.mfunction f=fun9(t,y)f=[y(2), 9.8*sin(y(1))]'; %f 向量必须为一列向量运行MATLAB 代码clear; close;[t,y]=ode45('fun9',[0,10],[15,0]);plot(t,y(:,1)); %画θ随时间变化图,y(:2)则表示'θ的值 xlabel('t'),ylabel('y1')结果见图6.7.2图6.7.2 数值解由图6.7.2可见，θ随时间t 周期变化。
练习4 (刚性方程组求解)求下面刚性微分方程的解⎪⎩⎪⎨⎧==-=--=1)0(,2)0(,100'99.9901.0'2122211y y y y y y y 使用dsolve 可知解析解为)100ex p(),100ex p()01.0ex p(21t y t t y -=-+-=下面求数值解. 建立M 文件fun10.m 如下%M 文件fun10.mfunction f=fun10(t,y)f=[-0.01*y(1)-99.99*y(2), -100*y(2)]';运行MATLAB 代码clear; close;[t,y]=ode45('fun10',[0,10],[2,1]);plot(t,y); text(1,1.1,'y1'); text(1,0.1,'y2'); xlabel('t'),ylabel('y')结果见图6.7.3图6.7.3 数值解图6.7.3给人的感觉似乎是1y 始终大于0.5.但由21,y y 的解析解可知,当∞→t 时,两个分量21,y y 均趋于0.2y 下降极快,0001.0)1.0(2<y ; 而1y 下降很慢,0183.0)400(2≈y (见下图6.7.4).若用clear; close;[t,y]=ode45('fun10',[0,400],[2,1]); tstep=length(t) %求计算总步数 minh=min(diff(t)) %最小步长 maxh=max(diff(t)) %最大步长结果为tstep =48261minh =5.0238e-004 maxh =0.0102可见计算太慢,t 需要48261步才能到达400.一方面,由于2y 下降太快,为了保证数值稳定性,步长h 须足够小;另一方面,由于1y 下降太慢,为了反映解的完整性,时间区间须足够长,这就造成计算量太大.这类方程称为刚性方程或病态方程.ode45不适用于病态方程,下面我们用ode15s 求解.clear; close;[t,y]=ode15s('fun10',[0,400],[2,1]);plot(t,y); text(100,0.5,'y1'); text(1,0.1,'y2'); xlabel('t'),ylabel('y') tstep=length(t) minh=min(diff(t)) maxh=max(diff(t)) 结果为tstep = 92minh =3.5777e-004 maxh =32.1282可见只需92步,最大步长为32,速度快了约500倍.函数图形见图6.7.4.图6.7.4 数值解练习5 (Lorenz 吸引子) 求常微分方程⎪⎪⎪⎩⎪⎪⎪⎨⎧+-=--=+-=xy z dtdz xz y x dtdyy x dt dx38281010 的数值解,初值取1)0()0()0(===z y x . 先建立M 文件Lorenzf.m 如下%M 文件Lorenzf.mfunction f=lorenzf(t,x)sig=10;bet=8/3;rho=28;f=[sig*(x(2)-x(1)),x(1).*(rho-x(3))-x(2),x(1).*x(2)-bet*x(3)]; f=f(:);运行MATLAB 代码clear; close;[t,y]=ode45('Lorenzf',[0,100],[1,1,1]); plot3(y(:,1),y(:,2),y(:,3)); xlabel x;ylabel y;zlabel z; 运行结果如图6.7.5.xyz图6.7.5 Lorenz 吸引子实验作业1．求下列微分方程的解析解(1) 一阶线性方程2'3=-y x y (2) 贝努利方程0'2=--y xy y(3) 高阶线性齐次方程02'3"'"=+--y y y y (4) 高阶线性非齐次方程x y y y sin 32'3"=+- 2．求方程3)0(',1)0(,'2")1(2===+y y xy y x的解析解和数值解,并进行比较3．分别用ode45和ode15s 求解Van-del-Pol 方程()⎪⎩⎪⎨⎧===---1)0',0)0(0)1(1000222x x x dt dxx dtx d 的数值解,并进行比较.4. （Rosseler 吸引子）用ode45数值求解方程⎪⎪⎪⎩⎪⎪⎪⎨⎧-+=+=+-=)()(c x z b dtdz ayx dtdyz y dt dx其中初值取8)0(,3)0(,2)0(===z y x ，参数7.5,2.0,2.0===c b a .阅读资料我国微分方程界的先辈申又枨教授申又枨，数学家、数学教育家。
从事函数论及微分方程的研究。
主要成就涉及复变函数的插值理论。
是在新中国建立微分方程学科研究的创始人之一。
申又枨，1901年6月13日生于山西高平鼓楼。