常微分方程的数值解
实验4 常微分方程的数值解
【实验目的】
1．掌握用MATLAB软件求微分方程初值问题数值解的方法；
2．通过实例用微分方程模型解决简化的实际问题；
3．了解欧拉方法和龙格-库塔方法的基本思想和计算公式，及稳定性等概念。
【实验内容】
题3
小型火箭初始重量为1400kg，其中包括1080kg燃料。火箭竖直向上发射时燃料燃烧率为18kg/s，由此产生32000N的推力，火箭引擎在燃料用尽时关闭。设火箭上升时空气阻力正比于速度的平方，比例系数为0.4kg/m，求引擎关闭瞬间火箭的高度、速度、加速度，及火箭到达最高点的时的高度和加速度，并画出高度、速度、加速度随时间变化的图形。
模型及其求解
火箭在上升的过程可分为两个阶段，在全过程中假设重力加速度始终保持不变，g=9.8m/s2。
在第一个过程中，火箭通过燃烧燃料产生向上的推力，同时它还受到自身重力（包括自重和该时刻剩余燃料的重量）以及与速度平方成正比的空气阻力的作用，根据牛顿第二定律，三个力的合力产生加速度，方向竖直向上。因此有如下二式：
a=dv/dt=(F-mg-0.4v2)/m=(32000-0.4v2)/(1400-18t)-9.8
dh/dt=v
又知初始时刻t=0，v=0，h=0。记x(1)=h，x(2)=v，根据MATLAB 可以求出0到60秒内火箭的速度、高度、加速度随时间的变化情况。程序如下：
function [ dx ] = rocket( t,x )
a=[(32000-0.4*x(2)^2)/(1400-18*t)]-9.8;
dx=[x(2);a];
end
ts=0:1:60;
x0=[0,0];
[t,x]=ode45(@rocket,ts,x0);
h=x(:,1);
v=x(:,2);
a=[(32000-0.4*(v.^2))./(1400-18*t)]-9.8; [t,h,v,a];
数据如下：
t h v a
0 0 0 13.06
1.00 6.57 13.19 13.30
2.00 26.44 26.58 1
3.45
3.00 59.76 40.06 13.50
4.00 106.57 53.54 13.43
5.00 16
6.79 66.89 13.26
6.00 240.27 80.02 12.99
7.00 326.72 92.83 12.61
8.00 425.79 105.22 12.15
9.00 536.99 117.11 11.62
10.00 659.80 128.43 11.02
11.00 793.63 139.14 10.38
12.00 937.85 149.18 9.71
13.00 1091.79 158.55 9.02
14.00 1254.71 167.23 8.33
15.00 1425.93 175.22 7.65
16.00 1604.83 182.55 6.99
17.00 1790.78 189.22 6.36
18.00 1983.13 195.27 5.76
19.00 2181.24 200.75 5.21
20.00 2384.47 205.70 4.69
21.00 2592.36 210.18 4.22
22.00 2804.52 214.19 3.79
23.00 3020.56 217.79 3.41
24.00 3240.08 221.01 3.07
25.00 3462.65 223.92 2.77
26.00 3687.88 226.56 2.50
27.00 3915.58 228.97 2.27
28.00 4145.60 231.14 2.06
29.00 4377.76 233.11 1.89
30.00 4611.86 234.91 1.74
31.00 4847.68 236.57 1.62
32.00 5085.02 238.14 1.51
33.00 5323.85 239.61 1.41
34.00 5564.11 240.99 1.33
35.00 5805.77 242.28 1.27
36.00 6048.72 243.50 1.21
37.00 6292.87 244.68 1.17
38.00 6538.11 245.83 1.13
39.00 6784.48 246.96 1.09
40.00 7031.96 248.05 1.07
41.00 7280.54 249.10 1.05
42.00 7530.19 250.12 1.03
43.00 7780.85 251.14 1.02
44.00 8032.49 252.15 1.00
45.00 8285.12 253.16 0.99
46.00 8538.75 254.15 0.98
47.00 8793.39 255.12 0.97
48.00 9049.01 256.07 0.97
49.00 9305.58 257.03 0.96
50.00 9563.08 257.99 0.95
51.00 9821.52 258.95 0.94
52.00 10080.93 259.90 0.93
53.00 10341.30 260.83 0.93
54.00 10602.62 261.75 0.94
55.00 10864.86 262.67 0.94
56.00 11127.98 263.61 0.93
57.00 11392.04 264.54 0.91
58.00 11657.03 265.46 0.91
59.00 11922.96 266.35 0.92
60.00 12189.78 267.26 0.92
因此，在引擎关闭的瞬间，火箭的速度为267.26m/s，高度为
12189.78m，加速度为0.92m/s2。
（2）在第二个阶段，火箭的重量保持不变，没有向上的推力，只收到重力和空气阻力。因此有如下关系式：
a=dv/dt=(-mg-0.4v2)/m=-0.4v2/320-9.8
dh/dt=v
假设在80秒之前达到最高点，以60秒时的速度、高度、加速度为初始值进行计算，程序如下：
function [ dx ] = rocket2( t,x )
dx=[x(2);-0.4*x(2)^2/320-9.8];
end
ts2=60:1:80;
x1=[12189.78,267.26];
[t2,x2]=ode45(@rocket2,ts2,x1);
h2=x2(:,1);
v2=x2(:,2);
a2=-0.4*v2.^2./320-9.8;
[t2,h2,v2,a2];
数据如下：
t2 h2 v2 a2
60.00 12189.78 267.26 -99.08
61.00 12416.32 192.70 -56.22
62.00 12584.73 147.43 -36.97
63.00 12715.67 116.11 -26.65
64.00 12819.59 92.79 -20.56
65.00 12902.81 74.29 -16.70
66.00 12969.22 58.96 -14.15
67.00 13021.42 45.73 -12.41
68.00 13061.17 33.95 -11.24
69.00 13089.63 23.12 -10.47
70.00 13107.61 12.91 -10.01
71.00 13115.55 3.02 -9.81
72.00 13113.66 -6.80 -9.86
73.00 13101.90 -16.78 -10.15
74.00 13079.96 -27.19 -10.72
75.00 13047.26 -38.34 -11.64
76.00 13002.90 -50.62 -13.00
77.00 12945.47 -64.56 -15.01
78.00 12872.96 -80.96 -17.99
79.00 12782.31 -101.08 -22.57
80.00 12668.88 -127.03 -29.97
可以看到：在第60秒瞬间，加速度发生了突变，从0.92m/s2突变为-99.08m/s2；而在第71秒至第72秒之间，速度从正变为负，即速度反向，说明在第71秒中的某个时刻速度为零，火箭达到了最高点。因此需要对这个时间段进行分析，并且找到速度减小到0的时刻和此时的高度。
以0.1为步长，在71s到72s中重新求解微分方程的数值解。
71.00 13115.16 2.98 -9.81
71.10 13115.41 2.00 -9.81
71.20 13115.56 1.02 -9.80
71.30 13115.61 0.04 -9.80
71.40 13115.57 -0.94 -9.80
71.50 13115.42 -1.92 -9.80
可见在t=71.3时，速度为0.04，可视为速度为零点，此时最大高度为13115.61，加速度-9.80。
综合（1），（2），可以绘出高度，速度和加速度随时间的变化曲线。
plot(t,h,t2,h2),xlabel('t/s'),ylabel('h/m'),title('高度随时间变化曲线');
plot(t,v,t2,v2),xlabel('t/s'),ylabel('v/(m/s)'),title('速度随时间变化曲线');
aa=[a',a2']'；tt=[t',t2']'；
plot(tt,aa),xlabel('t/s'),ylabel('a/(m/s2)'),title('加速度随时间的变化曲线');
题6
一只小船渡过宽为d 的河流，目标是起点A 正对着的另一岸B 点。已知河水流速v 1与船在静水中的速度v 2之比为k 。
（1）建立描述小船航线的数学模型，求其解析解；
（2）设d=100m ，v 1=1m/s ，v 2=2m/s ，用数值解法求渡河所需时间、任意时刻小船的位置及航行曲线，作图，并与解析解比较；
（3）若流速v 1=0，0.5，1.5，2m/s ，结果将如何。
模型及其求解
（1）.假设在航行过程中，人们不知道水流的速度，小船的方向始终指向目标B ，因此若以B 为原点，我们可以得到如下方程组：
解初值为（x, y ）=( 0, -d) 的常微分方程组，得到解析解为：
其中 c=–1/d=–0.01,故有
1v x dx
v dt v y dy dt ?=-??
?
?=-
??
1
1110
2
2
k k k
k
x c y
c
y
+--+
-
=1
111()(0.01)(0.01)
2
2
k
k k
k
x y y
y
+--=-
-+
-
事实上，若用matlab中计算微分方程的语句：
[x,y]=dsolve('Dx=v1-v2*x/sqrt(x^2+y^2)','Dy=-v2*y/sqrt(x^2+y^2)','x(0) =0','y(0)=-d','t');
则显示“Warning: Explicit solution could not be found.”即无法得到x,y关于t的分立解。
（2）d=100m,v1=1m/s,v2=2m/s。求数值解时，由解析解可以看出，此为刚性方程。为避免用ode45s求解时间过长。采用ode23s进行求解。假定100s 可以到达对岸。
ts=0:0.1:100;
x0=[0,-100];
[t,x]=ode23s(@boat,ts,x0);
[t,x];
plot(t,x);gtext('x');gtext('y');xlabel('时间t/s');ylabel('距离/m');title('x,y与时间t的关系');
可以得到数据如下（部分）：
t/s x/m y/m
10.00 8.93 -80.04
20.00 15.41 -60.36
30.00 18.87 -41.49
40.00 18.71 -24.31
50.00 14.48 -10.34
50.10 14.42 -10.23
66.50 0.19 -0.00
66.60 0.09 -0.00
66.70 0.00 0
可知在
t=66.7s时，船
到达对岸B点。
做x,y与时间t
的关系图：从曲
线上可以看出，0
到30s这段时间
内，y的增长几
乎呈线性关系，
即小船几乎研直
线匀速前进。
现在求解析解并将之与数值解对比：
function [ x] = jiexi(y,k)
x=-0.5*(-0.01)^k*y.^(k+1)+0.5*(-0.01)^(-k).*y.^(1-k);
end
y=-100:1:0;
x2=jiexi(y,0.5);
plot(x2,y,'ro',x(:,1),x(:,2));legend('解析解','数值解',-1);
从轨迹曲线中也可以看到，用数值解得到结果与解析解几乎重合，可信度很高。
（3）当改变v1的值时，解析式中的k值将发生变化，此时将对结果产生影响。利用MATLAB计算和绘图也可以发现，渡河时间及航行曲线都发生了变化。
v1=0时，k=0。说明河水静止不动，这种情况下，小船的总速度就是它在静水中的速度，于是沿着AB直线便可到达对岸，计算结果表示，t=50s时，小船到达B点。
这种情况下，t=53.33s时，小船到达B点，比起v1=1m/s时，小船在x方向上走过的距离缩短了大约一半，总路程缩短了许多。
这种情况下，t=114.28s，小船到达对岸，渡河时间明显长了许多。而且由数值解的疏密程度也可以看出，小船花费较少时间久到达x的最大位移处，但是却花了相当长的时间回到x=0的目的地，可见1.5m/s的河水流速给小船到达终点造成了巨大阻碍。
v1=2m/s时，k=1，得到的曲线如下：
这种情况下，小船无论如何都无法到达B点，只能到达B点下游50米处。从解析式中也可以看到，k=1时，有x(y)=-0.005y2+50。曲线呈开口朝左的抛物线状。从速度合成的三角形上来看，由于v1和v2长度相等，v1的方向也已确定，它们的合速度的方向与v1的夹角不可能大于90°，也就是说在x的分位移始终在增加，不可能减小。即使v2沿着水流逆向，也只能使合速度为零，此时正好是小船停在B点下游50米处的情况。类似的问题在高中物理出现过。
综合上述曲线，有下图：
仔细观察解析式可以发现影响船过河轨迹仅是k值，即船与河水的相对速度，我们不妨作此假设。
如果我们用v1=2m/s，v2=4m/s与前面的结果做下对比（我们仅做数值解的对比，因为解析解相同）。结果证明当v1=2m/s，v2=4m/s 时船的航行轨迹与v1=1m/s，v2=2m/s 航行轨迹相同，但时间缩短为33.33s。而且继续实验当
v1=4m/s，v2=8m/s时，船的航行轨迹也与前两种情况相同，但过河时间为16.66s，说明过河时间与船速成倒数关系，这是从航行轨迹相同、路程相等可以很自然导出的关系。
在实际问题中，人们通常会对水的流速做初步判断，以调整船行驶的方向，而不是单纯的每个时刻都将方向对准目的地。同时由于水的流速也不是始终不变的，在不同的流域和不同的时刻流速都可能不同，本题中只是采取了最为简单的数学模型进行近似计算。
题9两种群相互竞争的模型如下：
X’(t)=r1*x*(1-x/n1-s1*y/n2);
Y’(t)=r2*y*(1-s2*x/n1-y/n2);
其中X(t)，Y(t)分别为甲乙两种种群的数量；r1,r2为他们的固有增长率；n1，n2为他们的最大容量。s1的含义是：对于供养甲的资源而言，单位数量乙（相对n2）的消耗为单位数量的甲（相对n1）消耗的s1倍，对s2可做相应的解释。
该模型无解析解，使用数值的解法研究问题：
（1）设r1=r2=1，n1=n2=100，s1=0.5，s2=2，初值x0=y0=10，计算X(t)，Y(t)，画出他们的图形及相图（x,y），说明时间t充分大以后X(t)，Y(t)的变化趋势。
（2）改变r1，r2，n1，n2，x0，y0，但保持s1，s2不变（或保持s1<1，s2>1），计算并分析所得结果；若s1=1.5，s2=0.7，再分析结果。由此你能得到什么结论，请用各参数生态学上的含义做出解释。
（3）实验s1=0.8，s2=0.7和s1=1.5，s2=1.7时会有什么结果，并作出解释。
模型及其求解
（1）编写程序如下：
function [ dx ] = compt( t,x )
r1=1;r2=1;n1=100;n2=100;s1=0.5;s2=2;
dx=[r1*x(1)*(1-x(1)/n1-s1*x(2)/n2);r2*x(2)*(1-s2*x(1)/n1-x(2)/n2)];
end
ts=0:0.5:20;
x0=[10,10];
[t,x]=ode23s(@compt,ts,x0);
[t,x];
plot(t,x);title('t充分大后x(t),y(t)的变化趋势');xlabel('时间t');ylabel('种群数量');gtext('x(t)');gtext('y(t)');
plot(x(:,1),x(:,2));xlabel('x');ylabel('y');title('相图x,y');
设定t从0到20，得出t,x(t),y(t)的数值关系
t X(t) Y(t)
0 10.0000 10.0000
0.5000 15.0556 13.7363
1.0000 21.8041 17.4479
1.5000 30.1504 20.2070
2.0000 39.6707 21.1843
2.5000 49.6905 20.1332
3.0000 59.4906 17.4859
3.5000 68.4731 1
4.0317
4.0000 76.2336 10.5329
4.5000 82.5937 7.4878
5.0000 87.5743 5.0991
5.5000 91.3224 3.3592
6.0000 94.0516 2.1579
6.5000 95.9832 1.3614
7.0000 97.3207 0.8479
7.5000 98.2306 0.5233
8.0000 98.8409 0.3209
8.5000 99.2457 0.1960
9.0000 99.5118 0.1194
9.5000 99.6855 0.0726
10.0000 99.7982 0.0441
10.5000 99.8709 0.0267
11.0000 99.9177 0.0162
11.5000 99.9477 0.0098
12.0000 99.9668 0.0060
12.5000 99.9790 0.0036
13.0000 99.9867 0.0022
13.5000 99.9916 0.0013
14.0000 99.9947 0.0008
14.5000 99.9967 0.0005
15.0000 99.9979 0.0003
15.5000 99.9987 0.0002
16.0000 99.9992 0.0001
16.5000 99.9995 0.0001
17.0000 99.9997 0.0000
17.5000 99.9998 0.0000
18.0000 99.9999 0.0000
18.5000 99.9999 0.0000
19.0000 100.0000 0.0000 19.5000 100.0000 0.0000
20.0000 100.0000 0.0000
由统计数据可以看出，t=17时，乙种群的数量已经为0，之后甲种群的数量达到饱和。由图线更能很好的分析出甲乙两个种群的变化趋势。刚开始时两种群的数量都在增长，但环境的容纳量有限，甲乙两种群不可避免的竞争，于是甲的数量持续增长而乙数量在达到一个峰值后逐渐减少，两种群数量相差悬殊。最后乙灭绝而甲繁荣。
究其原因，s1 （2）探究r1，r2，n1，n2，x0，y0对种群数量的影响。除探究的因素以外，要让其他因素保持（1）中数值不变。 A. 首先实验自然增长率r。对于r1>r2，显然甲增加的会更快，结果与（1）相同。而对于r2>r1，不妨设r1=1,r2=4,作图。 可见开始的时候，凭借着高自然增长率，乙种群数量超过了甲。但之后甲种群数量开始稳步上升，乙种群数量则不可避免的下降，最后甲种群再次繁荣，乙种群灭绝。且达到最后稳定的时间比（1）要短，也许是食物消耗过快的结果。 由此可见，乙自然增长率r的提高并不能让乙种群避免灭绝的命运。 B. 其次探究最大容量n的影响。显而易见，如果n1>n2，则趋势与（1）相同。仅仅是甲的最大数量更大了。 对于n2>n1，不妨设n2=300,n1=100。作图： 可见乙种群最大容量的增大，也不能避免他灭绝的命运。从生态学上来看，这样的结果符合常理。 C. 之后探究种群数量初值。对于x0>y0，显然趋势和最后结果与（1）相同，甲的增长会更加迅速。对于x0 可见，乙种群初始数量的增大，虽然在开始时有过一个高峰，但是并不能避免最后灭绝的命运。从生态学上解释，当一个种群生存能力很弱时，即使初始数量很多，最后也会灭绝。 下面探究s1，s2。设s1=1.5。s2=0.7。其他因素与（1）相同。 可见结果发生了逆转：甲种群最终走向了衰落，几乎灭绝，而乙种群走向了繁荣。 由此看出，能否消耗更多供养自己和供养对手的食物和资源，是一个种群能否繁衍生存的关键因素。如果自己的食物和对手的食物都不能有效地“为我所用”，那么一个种群就难逃灭绝的命运。 数据如下： t X(t) Y(t) 0 10.0000 10.0000 0.5000 14.1649 14.8747 1.0000 18.8098 21.1830 1.5000 23.1789 28.6762 2.0000 26.4025 36.8127 2.5000 27.9436 44.9840 3.0000 27.7660 52.6886 3.5000 26.2263 59.6679 常微分方程数值解法 【作用】 微分方程建模是数学建模的重要方法，因为许多实际问题的数学描述将导致求解微 分方程的定解问题。把形形色色的实际问题化成微分方程的定解问题，大体上可以按以 下几步： 1. 根据实际要求确定要研究的量(自变量、未知函数、必要的参数等)并确定坐标系。 2. 找出这些量所满足的基本规律(物理的、几何的、化学的或生物学的等等)。 3. 运用这些规律列出方程和定解条件。 基本模型 1.发射卫星为什么用三级火箭 2.人口模型 3.战争模型 4.放射性废料的处理 通常需要求出方程的解来说明实际现象，并加以检验。如果能得到解析形式的解固然是便于分析和应用的，但是我们知道，只有线性常系数微分方程，并且自由项是某些特殊类型的函数时，才可以得到这样的解，而绝大多数变系数方程、非线性方程都是所谓“解不出来”的于是对于用微分方程解决实际问题来说，数值解法就是一个十分重要的手段。 1.改进Euler 法： 2.龙格—库塔（Runge—Kutta）方法： 【源程序】 1.改进Euler 法： function [x,y]=eulerpro(fun,x0,x1,y0,n);%fun为函数，(x0，x1)为x区间， y0为初始值，n为子区间个数 if nargin<5,n=50;end h=(x1-x0)/n; x(1)=x0;y(1)=y0; for i=1:n x(i+1)=x(i)+h; y1=y(i)+h*feval(fun,x(i),y(i)); y2=y(i)+h*feval(fun,x(i+1),y1); y(i+1)=(y1+y2)/2; end 调用command窗口 f=inline('-2*y+2*x^2+2*x') [x,y]=eulerpro(f,0,0.5,1,10) 求解函数y'=?2y+2x2+2x ，(0 ≤x ≤0.5)，y(0) = 1 2.龙格—库塔（Runge—Kutta）方法： [t,y]=solver('F',tspan，y0) 这里solver为ode45，ode23，ode113，输入参数 F 是用M文件定义的微分方程y'= f (x, y) 右端的函数。tspan=[t0,tfinal]是求解区间，y0是初值。 注：ode45和ode23变步长的，采用Runge-Kutta算法。 第四章常微分方程数值解 [课时安排]6学时 [教学课型]理论课 [教学目的和要求] 了解常微分方程初值问题数值解法的一些基本概念，如单步法和多步法，显式和隐式，方法的阶数，整体截断误差和局部截断误差的区别和关系等；掌握一阶常微分方程初值问题的一些常用的数值计算方法，例如欧拉（Euler）方法、改进的欧拉方法、龙贝－库塔（Runge-Kutta）方法、阿达姆斯（Adams）方法等，要注意各方法的特点及有关的理论分析；掌握构造常微分方程数值解的数值积分的构造方法和泰勒展开的构造方法的基本思想，并能具体应用它们导出一些常用的数值计算公式及评估截断误差；熟练掌握龙格－库塔（Ｒ－Ｋ）方法的基本思想，公式的推导，Ｒ－Ｋ公式中系数的确定，特别是能应用“标准四阶Ｒ－Ｋ公式”解题；掌握数值方法的收敛性和稳定性的概念，并能确定给定方法的绝对稳定性区域。[教学重点与难点] 重点：欧拉方法，改进的欧拉方法，龙贝－库塔方法。 难点：R—K方法，预估-校正公式。 [教学内容与过程] 4.1 引言 本章讨论常微分方程初值问题 (4.1.1) 的数值解法，这也是科学与工程计算经常遇到的问题，由于只有很特殊的方程能用解析方法求解，而用计算机求解常微分方程的初值问题都要采用数值方法.通常我们假定(4.1.1)中 f(x,y)对y满足Lipschitz条件，即存在常数L＞0，使对，有 (4.1.2) 则初值问题(4.1.1)的解存在唯一. 假定(4.1.1)的精确解为，求它的数值解就是要在区间上的一组离散点 上求的近似.通常取 ，h称为步长，求(4.1.1)的数值解是按节点的顺序逐步 推进求得.首先，要对方程做离散逼近，求出数值解的公式，再研究公式的局部截 第七章 常微分方程初值问题的数值解法 --------学习小结 一、本章学习体会 通过本章的学习，我了解了常微分方程初值问题的计算方法，对于解决那些很难求解出解析表达式的，甚至有解析表达式但是解不出具体的值的常微分方程非常有用。在这一章里求解常微分方程的基本思想是将初值问题进行离散化，然后进行迭代求解。在这里将初值问题离散化的方法有三种，分别是差商代替导数的方法、Taylor 级数法和数值积分法。常微分方程初值问题的数值解法的分类有显示方法和隐式方法，或者可以分为单步法和多步法。在这里单步法是指计算第n+1个y 的值时,只用到前一步的值，而多步法则是指计算第n+1个y 的值时，用到了前几步的值。通过对本章的学习，已经能熟练掌握如何用Taylor 级数法去求解单步法中各方法的公式和截断误差，但是对线性多步法的求解理解不怎么透切，特别是计算过程较复杂的推理。 在本章的学习过程中还遇到不少问题，比如本章知识点多，公式多，在做题时容易混淆，其次对几种R-K 公式的理解不够透彻，处理一个实际问题时，不知道选取哪一种公式，通过课本里面几种方法的计算比较得知其误差并不一样，，这个还需要自己在往后的实际应用中多多实践留意并总结。 二、本章知识梳理 常微分方程初值问题的数值解法一般概念 步长h ，取节点0,(0,1,...,)n t t nh n M =+=，且M t T ≤，则初值问题000 '(,),()y f t y t t T y t y =≤≤?? =?的数值解法的一般形式是 1(,,,...,,)0,(0,1,...,)n n n n k F t y y y h n M k ++==- 第九章 常微分方程的数值解法 在自然科学的许多领域中，都会遇到常微分方程的求解问题。然而，我们知道，只有少数十分简单的微分方程能够用初等方法求得它们的解，多数情形只能利用近似方法求解。在常微分方程课中已经讲过的级数解法，逐步逼近法等就是近似解法。这些方法可以给出解的近似表达式，通常称为近似解析方法。还有一类近似方法称为数值方法，它可以给出解在一些离散点上的近似值。利用计算机解微分方程主要使用数值方法。 我们考虑一阶常微分方程初值问题 ?????==0 0)() ,(y x y y x f dx dy 在区间[a, b]上的解，其中f (x, y )为x, y 的已知函数，y 0为给定的初始值，将上述问题 的精确解记为y (x )。数值方法的基本思想是：在解的存在区间上取n + 1个节点 b x x x x a n =<<<<= 210 这里差i i i x x h -=+1，i = 0,1, …, n 称为由x i 到x i +1的步长。这些h i 可以不相等，但一 般取成相等的，这时n a b h -=。在这些节点上采用离散化方法，（通常用数值积分、微分。泰勒展开等）将上述初值问题化成关于离散变量的相应问题。把这个相应问题的解y n 作为y (x n )的近似值。这样求得的y n 就是上述初值问题在节点x n 上的数值解。一般说来，不同的离散化导致不同的方法。 §1 欧拉法与改进欧拉法 1.欧拉法 1．对常微分方程初始问题 (9.2) )((9.1) ),(00???? ?==y x y y x f dx dy 用数值方法求解时，我们总是认为(9.1)、(9.2)的解存在且唯一。 欧拉法是解初值问题的最简单的数值方法。从（9.2）式由于y (x 0) = y 0已给定，因而可以算出 ),()('000y x f x y = 设x 1 = h 充分小，则近似地有： ),()(') ()(00001y x f x y h x y x y =≈- (9.3) i.常微分方程初值问题数值解法 常微分方程初值问题的真解可以看成是从给定初始点出发的一条连续曲线。差分法是常微分方程初值问题的主要数值解法，其目的是得到若干个离散点来逼近这条解曲线。有两个基本途径。一个是用离散点上的差商近似替代微商。另一个是先对微分方程积分得到积分方程，再利用离散点作数值积分。 i.1 常微分方程差分法 考虑常微分方程初值问题：求函数()u t 满足 (,), 0du f t u t T dt =<≤ （i.1a ） 0(0)u u = (i.1b) 其中(,)f t u 是定义在区域G : 0t T ≤≤, u <∞上的连续函数，0u 和T 是给定的常数。我们假设(,)f t u 对u 满足Lipschitz 条件，即存在常数L 使得 121212(,)(,), [0,]; ,(,)f t u f t u L u u t T u u -≤-?∈∈-∞∞ （i.2） 这一条件保证了（i.1）的解是适定的，即存在，唯一，而且连续依赖于初值0u 。 通常情况下，（i.1）的精确解不可能用简单的解析表达式给出，只能求近似解。本章讨论常微分方程最常用的近似数值解法-差分方法。先来讨论最简单的Euler 法。为此，首先将求解区域[0,]T 离散化为若干个离散点： 0110N N t t t t T -=<< <<= （i.3） 其中n t hn =，0h >称为步长。 在微积分课程中我们熟知，微商（即导数）是差商的极限。反过来，差商就是微商的近似。在0t t =处，在（i.1a ）中用向前差商 10()()u t u t h -代替微商du dt ，便得 10000()()(,())u t u t hf t u t ε=++ 如果忽略误差项0ε，再换个记号，用i u 代替()i u t 便得到 1000(,)u u hf t u -= 一般地，我们有 1Euler (,), 0,1, ,1n n n n u u hf t u n N +=+=-方法： （i.4） 从(i.1b) 给出的初始值0u 出发，由上式可以依次算出1,,N t t 上的差分解1,,N u u 。 利用MATLAB求解常微分方程数值解 目录 1. 内容简介 (1) 2. Euler Method(欧拉法)求解 (1) 2.1. 显式Euler法和隐式Euler法 (2) 2.2. 梯形公式和改进Euler法 (3) 2.3. Euler法实用性 (4) 3. Runge-Kutta Method(龙格库塔法)求解 (5) 3.1. Runge-Kutta基本原理 (5) 3.2. MATLAB中使用Runge-Kutta法的函数 (7) 4. 使用MATLAB求解常微分方程 (7) 4.1. 使用ode45函数求解非刚性常微分方程 (8) 4.2. 刚性常微分方程 (9) 5. 总结 (9) 参考文献 (11) 附录 (12) 1. 显式Euler法数值求解 (12) 2. 改进Euler法数值求解 (12) 3. 四阶四级Runge-Kutta法数值求解 (13) 4.使用ode45求解 (14) 1.内容简介 把《高等工程数学》看了一遍，增加对数学内容的了解，对其中数值解法比较感兴趣，这大概是因为在其它各方面的学习和研究中经常会遇到数值解法的问题。理解模型然后列出微分方程，却对着方程无从下手，无法得出精确结果实在是让人难受的一件事情。 实际问题中更多遇到的是利用数值法求解偏微分方程问题，但考虑到先从常微分方程下手更为简单有效率，所以本文只研究常微分方程的数值解法。把一个工程实际问题弄出精确结果远比弄清楚各种细枝末节更有意思，因此文章中不追求非常严格地证明，而是偏向如何利用工具实际求解出常微分方程的数值解，力求将课程上所学的知识真正地运用到实际方程的求解中去，在以后遇到微分方程的时候能够熟练运用MATLAB得到能够在工程上运用的结果。 文中求解过程中用到MATLAB进行数值求解，主要目的是弄清楚各个函数本质上是如何对常微分方程进行求解的，对各种方法进行MATLAB编程求解，并将求得的数值解与精确解对比，其中源程序在附录中。最后考察MATLAB中各个函数的适用范围，当遇到实际工程问题时能够正确地得到问题的数值解。 2.Euler Method(欧拉法)求解 Euler法求解常微分方程主要包括3种形式，即显式Euler法、隐式Euler法、梯形公式法，本节内容分别介绍这3种方法的具体内容，并在最后对3种方法精度进行对比，讨论Euler法的实用性。 本节考虑实际初值问题 使用解析法，对方程两边同乘以得到下式 实验4 常微分方程的数值解 【实验目的】 1．掌握用MATLAB软件求微分方程初值问题数值解的方法； 2．通过实例用微分方程模型解决简化的实际问题； 3．了解欧拉方法和龙格-库塔方法的基本思想和计算公式，及稳定性等概念。 【实验内容】 题3 小型火箭初始重量为1400kg，其中包括1080kg燃料。火箭竖直向上发射时燃料燃烧率为18kg/s，由此产生32000N的推力，火箭引擎在燃料用尽时关闭。设火箭上升时空气阻力正比于速度的平方，比例系数为m，求引擎关闭瞬间火箭的高度、速度、加速度，及火箭到达最高点的时的高度和加速度，并画出高度、速度、加速度随时间变化的图形。 模型及其求解 火箭在上升的过程可分为两个阶段，在全过程中假设重力加速度始终保持不变，g=s2。 在第一个过程中，火箭通过燃烧燃料产生向上的推力，同时它还受到自身重力（包括自重和该时刻剩余燃料的重量）以及与速度平方成正比的空气阻力的作用，根据牛顿第二定律，三个力的合力产生加速度，方向竖直向上。因此有如下二式： a=dv/dt=/m=/(1400-18t) dh/dt=v 又知初始时刻t=0，v=0，h=0。记x(1)=h，x(2)=v，根据MATLAB 可以求出0到60秒内火箭的速度、高度、加速度随时间的变化情况。程序如下： function [ dx ] = rocket( t,x ) a=[*x(2)^2)/(1400-18*t)]; dx=[x(2);a]; end ts=0:1:60; x0=[0,0]; [t,x]=ode45(@rocket,ts,x0); h=x(:,1); v=x(:,2); a=[*(v.^2))./(1400-18*t)]; [t,h,v,a]; 数据如下： t h v a 000 第七章 常微分方程数值解法 常微分方程中只有一些典型方程能求出初等解(用初等函数表示的解)，大部分的方程是求不出初等解的。另外，有些初值问题虽然有初等解，但由于形式太复杂不便于应用。因此，有必要探讨常微分方程初值问题的数值解法。本章主要介绍一阶常微分方程初值问题的欧拉法、龙格-库塔法、阿达姆斯方法，在此基础上推出一阶微分方程组与高阶方程初值问题的 数值解法；此外，还将简要介绍求解二阶常微分方程值问题的差分方法、试射法。 第一节 欧拉法 求解常微分方程初值问题 ?????==0 0)() ,(y x y y x f dx dy （1） 的数值解，就是寻求准确解)(x y 在一系列离散节点 <<<< 常微分方程初值问题数值解法 朱欲辉 （浙江海洋学院数理信息学院, 浙江舟山316004） [摘要]：在常微分方程的课程中讨论的都是对一些典型方程求解析解的方法.然而在生产实 际和科学研究中所遇到的问题往往很复杂, 在很多情况下都不可能给出解的解析表达式. 本篇文章详细介绍了常微分方程初值问题的一些数值方法, 导出了若干种数值方法, 如Euler法、改进的Euler法、Runge－Kutta法以及线性多步法中的Adams显隐式公式和预测校正 公式, 并且对其稳定性及收敛性作了理论分析. 最后给出了数值例子, 分别用不同的方法计算出近似解, 从得出的结果对比各种方法的优缺点. [关键词]：常微分方程;初值问题; 数值方法; 收敛性; 稳定性; 误差估计 Numerical Method for Initial-Value Problems Zhu Yuhui (School of Mathematics, Physics, and Information Science, Zhejiang Ocean University, Zhoushan, Zhejiang 316004) [Abstract]:In the course about ordinary differential equations, the methods for analytic solutions of some typical equations are often discussed. However, in scientific research, the problems are very complex and the analytic solutions about these problems can’t be e xpressed explicitly. In this paper, some numerical methods for the initial-value problems are introduced. these methods include Euler method, improved Euler method, Runge-Kutta method and some linear multistep method (e.g. Adams formula and predicted-corrected formula). The stability and convergence about the methods are presented. Some numerical examples are give to demonstrate the effectiveness and accuracy of theoretical analysis. [Keywords]:Ordinary differential equation; Initial-value problem; Numerical method; Convergence; Stability;Error estimate 第八章 常微分方程的数值解法 一．内容要点 考虑一阶常微分方程初值问题:?????==0 0)() ,(y x y y x f dx dy 微分方程的数值解：设微分方程的解y (x )的存在区间是[a,b ]，在[a,b ]内取一系列节 点a= x 0< x 1<…< x n =b ,其中h k =x k+1-x k ；(一般采用等距节点，h=(b-a)/n 称为步长)。在每个节点x k 求解函数y(x)的近似值:y k ≈y(x k )，这样y 0 , y 1 ,...,y n 称为微分方程的数值解。 用数值方法，求得f(x k )的近似值y k ，再用插值或拟合方法就求得y(x)的近似函数。 (一)常微分方程处置问题解得存在唯一性定理 对于常微分方程初值问题:?????==0 0)() ,(y x y y x f dx dy 如果： (1) 在B y y A x x 00≤-≤≤,的矩形内),(y x f 是一个二元连续函数。 (2) ),(y x f 对于y 满足利普希茨条件，即 2121y y L y x f y x f -≤-),(),(则在C x x 0≤≤上方程?????==0 0)() ,(y x y y x f dx dy 的解存在且唯一，这里C=min((A-x 0),x 0+B/L),L 是利普希茨常数。 定义：任何一个一步方法可以写为),,(h y x h y y k k k 1k Φ+=+，其中),,(h y x k k Φ称为算法的增量函数。 收敛性定理：若一步方法满足： (1)是p 解的. (2) 增量函数),,(h y x k k Φ对于y 满足利普希茨条件. (3) 初始值y 0是精确的。则),()()(p h O x y kh y =-kh =x -x 0,也就是有 0x y y lim k x x kh 0h 0 =--=→)( (一)、主要算法 1．局部截断误差 局部截断误差：当y(x k )是精确解时，由y(x k )按照数值方法计算出来的1~ +k y 的误差y (x k+1)- 1~ +k y 称为局部截断误差。 注意：y k+1和1~ +k y 的区别。因而局部截断误差与误差e k +1=y (x k +1) -y k +1不同。 如果局部截断误差是O (h p+1),我们就说该数值方法具有p 阶精度。 i.常微分方程初值问题数值解法 i.1 常微分方程差分法 考虑常微分方程初值问题：求函数()u t 满足 (,), 0du f t u t T dt =<≤ （i.1a ） 0(0)u u = (i.1b) 其中(,)f t u 是定义在区域G : 0t T ≤≤, u <∞上的函数，0u 和T 是给定的常数。我们假设(,)f t u 对u 满足Lipschitz 条件，即存在常数L 使得 121212(,)(,), [0,]; ,(,)f t u f t u L u u t T u u -≤-?∈∈-∞∞ （i.2） 这一条件保证了（i.1）的解是适定的，即存在，唯一，而且连续依赖于初值0u 。 通常情况下，（i.1）的精确解不可能用简单的解析表达式给出，只能求近似解。本章讨论常微分方程最常用的近似数值解法--差分方法。先来讨论最简单的Euler 法。为此，首先将求解区域[0,]T 离散化为若干个离散点： 0110N N t t t t T -=<<<<=L （i.3） 其中n t hn =，0h >称为步长。 在微积分课程中我们熟知，微商（即导数）是差商的极限。反过来，差商就是微商的近似。在0t t =处，在（i.1a ）中用向前差商 10()()u t u t h -代替微商du dt ，便得 10000()()(,())u t u t hf t u t ε=++ 如果忽略误差项0ε，再换个记号，用i u 代替()i u t 便得到 1000(,)u u hf t u -= 一般地，我们有 1Euler (,), 0,1,,1n n n n u u hf t u n N +=+=-L 方法： （i.4） 从(i.1b) 给出的初始值0u 出发，由上式可以依次算出1,,N t t L 上的差分解1,,N u u L 。 第一章绪论 1.1 引言 常微分方程是现代数学的一个重要分支，是人们解决各种实际问题的有效工具。微分方程的理论和方法从17世纪末开始发展起来，很快成了研究自然现象的强有力工具，在17到18世纪，在力学、天文、科学技术、物理中，就已借助微分方程取得了巨大的成就。1864年Leverrer根据这个方程预见了海王星的存在，并确定出海王星在天空中的位置。现在，常微分方程在许多方面获得了日新月异的应用。这些应用也为常微分方程的进一步发展提供了新的问题，促使人们对微分方程进行更深入的研究，以便适应科学技术飞速发展的需要。 研究常微分方程常用数值解是数学工作者的一项基本的且重要的工作。在国内外众多数学家的不懈努力，使此学科基本上形成了一套完美的体系。微分方程的首要问题是如何求一个给定方程的通解或特解。到目前为止，人们已经对许多微分方程得出了求解的一般方法。由于在生产实际和科学研究中所遇到的微分方程问题比较复杂，使这些问题的解即使能求出解析表达式，也往往因计算量太大而难于求出，而对于一些典型的微分方程则可以运用基本方法求出其解析解，并可以根据初值问题的条件把其中的任意常数确定下来。 由于求通解存在许多困难，人们就开始研究带某种定解条件的特解。首先是Cauchy对微分方程初始解的存在惟一性进行了研究。目前解的存在惟一性、延拓性、大范围的存在性以及解对初始解和参数的延续性和可微性等理论问题都已发展成熟。与此同时，人们开始采取各种近似方法来求微分方程的特解，例如求微分方程数值解的Euler折线法、Runge-Kutta法等，可以求得若干个点上微分方程的近似解。最后，由于当代高科技的发展为数学的广泛应用和深入研究提供了更好的手段。用计算机结合Matlab软件求方程的精确解、近似解，对解的性态进行图示和定性、稳定性研究都十分方便有效。 本章先介绍常微分的一般概念、导出微分方程的一些典型例子及求解微分方程的思路分析。从而得到常微分方程的常用数值解法。 第八章 常微分方程数值解法 考核知识点： 欧拉法，改进欧拉法，龙格-库塔法，单步法的收敛性与稳定性。 考核要求： 1． 解欧拉法，改进欧拉法的基本思想；熟练掌握用欧拉法，改进欧拉法、求微 分方程近似解的方法。 2． 了解龙格-库塔法的基本思想；掌握用龙格-库塔法求微分方程近似解的方 法。 3． 了解单步法的收敛性、稳定性与绝对稳定性。 例1 用欧拉法，预估——校正法求一阶微分方程初值问题 ? ??=-='1)0(y y x y ，在0=x （0，1）0.2近似解 解 （1）用1.0=h 欧拉法计算公式 n n n n n n x y y x y y 1.09.0)(1.01+=-+=+，1.0=n 计算得 9.01=y 82.01.01.09.09.02=?+?=y （2）用预估——校正法计算公式 1,0)(05.01.09.0)0(111)0(1=???-+-+=+=++++n y x y x y y x y y n n n n n n n n n 计算得 91.01=y ，83805.02=y 例2 已知一阶初值问题 ???=-='1 )0(5y y y 求使欧拉法绝对稳定的步长n 值。 解 由欧拉法公式 n n n n y h y h y y )51(51-=-=+ n n y h y ~)51(~1-=+ 相减得01)51()51(e h e h e n n n -==-=-Λ 当 151≤-h 时，4.00≤ Ch5. 常微分方程数值解法 §1. 引言 1. 问题的提出 假设一阶微分方程初值问题???=='00 )() ,(y x y y x f y 中的(),f x y 关于y 满足Lipschitz 条件，即存在常数L ，使得()()1212,,f x y f x y L y y -≤-，则由常微分方程理论知，初值问题有唯一解。 除了一些特殊类型的方程外，许多微分方程都没有解析解。 2. 数值解法的基本思想——离散化 计算解)(x y 在离散点 ,,,,10n x x x 上值)(i x y 的近似值i y ，ih x x i +=0。 3. 几个基本概念 (1) 单步法与多步法 若计算1i y +时只用到i y ，则称这种方法为单步法，如()1,i i i i y y hf x y +=+；若计算1i y +时需用到()1,, ,1i i i k y y y k --≥，则称这种方法为多步法。 (2) 显式与隐式 若1i y +可以直接用1,,,i i i k y y y --表示，则称此计算公式为显式，否则称之为 隐式。 §2. Euler 方法 1. Euler 公式 将),(y x f y ='在[]1,+n n x x 上积分，?+=-+1))(,()()(1n n x x n n dx x y x f x y x y ， 得? +=≈-+1))(,(1n n x x n n I dx x y x f y y ，用数值积分法求I 。 (1) ()n n y x hf I ,=，得()n n n n y x hf y y ,1+=+。 Euler 公式 (2) ()11,++=n n y x hf I ，得()111,++++=n n n n y x hf y y 。 后退的Euler 公式 (3) ()()[]11,,2 +++=n n n n y x f y x f h I ， 得()()[]111 ,,2 +++++=n n n n n n y x f y x f h y y 。 梯形公式（隐式） 数值分析--第9章常微分方程数值解 第九章 常微分方程数值解法 许多实际问题的数学模型是微分方程或微分方程的定解问题。如物体运动、电路振荡、化学反映及生物群体的变化等。常微分方程可分为线性、非线性、高阶方程与方程组等类；线性方程包含于非线性类中，高阶方程可化为一阶方程组。若方程组中的所有未知量视作一个向量，则方程组可写成向量形式的单个方程。因此研究一阶微分方程的初值问题 ?????=≤≤=0)(),(y a y b x a y x f dx dy ， (9-1) 的数值解法具有典型性。 常微分方程的解能用初等函数、特殊函数或它们的级数与积分表达的很少。用解析方法只能求出线性常系数等特殊类型的方程的解。对非线性方程来说，解析方法一般是无能为力的，即使某些解具有解析表达式，这个表达式也可能非常复杂而不便计算。因此研究微分方程的数值解法是非常必要的。 只有保证问题(9-1)的解存在唯一的前提下，研究其数值解法或者说寻求其数值解才有意义。由常微分方程的理论知，如果(9-1)中的),(y x f 满足条件 （1）),(y x f 在区域} ),({+∞<<∞-≤≤=y b x a y x D ，上连续； （2）),(y x f 在D 上关于y 满足Lipschitz 条件，即存在常数L ，使得 y y L y x f y x f -≤-),(),( 则初值问题(9-1)在区间],[b a 上存在惟一的连续解)(x y y =。在 下面的讨论中，我们总假定方程满足以上两个条件。 所谓数值解法，就是求问题(9-1)的解)(x y y =在若干点 b x x x x a N =<<<<= 210 处的近似值),,2,1(N n y n =的方法。),,2,1(N n y n =称为问题(9-1)的 数值解，n n x x h -=+1称为由n x 到1+n x 的步长。今后如无特别说明， 我们总假定步长为常量。 建立数值解法，首先要将微分方程离散化，一般采用以下几种方法： (1) 用差商近似导数 在问题(9-1)中，若用向前差商h x y x y n n )() (1-+代替)(n x y '，则得 )1,,1,0( ))(,()()(1-=≈-+N n x y x f h x y x y n n n n n )(n x y 用其近似值n y 代替，所得结果作为)(1+n x y 的近似值，记为1+n y ，则有 1(,) (0,1,,1)n n n n y y hf x y n N +=+=- 这样，问题(9-1)的近似解可通过求解下述问题 100 (,) (0,1,,1)()n n n n y y hf x y n N y y x +=+=-??=? (9-2) 得到，按式(9-2)由初值0y 经过N 步迭代，可逐次算出N y y y ,,21。此方程称为差分方程。 需要说明的是，用不同的差商近似导数，将得到不同的计算公式。 (2) 用数值积分法 将问题(9-1)中的微分方程在区间],[1+n n x x 上两边积分，可得 )1,,1,0( ))(,()()(1 1-==-?++N n dx x y x f x y x y n n x x n n (9-3) 用1+n y ,n y 分别代替)(1+n x y ,)(n x y ，若对右端积分采用取左端点的矩形公式，即 2.4 常微分方程数值解 函数 ode45、ode23、ode113、ode15s 、ode23s 、ode23t 、ode23tb 功能 常微分方程（ODE ）组初值问题的数值解 参数说明： solver 为命令ode45、ode23,ode113,ode15s,ode23s,ode23t,ode23tb 之一。 Odefun 为显式常微分方程y’=f(t,y)，或为包含一混合矩阵的方程M(t,y)*y’=f(t,y)。命令 ode23只能求解常数混合矩阵的问题；命令ode23t 与ode15s 可以求解奇异矩阵的问题。 Tspan 积分区间（即求解区间）的向量tspan=[t0,tf]。要获得问题在其他指定时间点 t0,t1,t2,…上的解，则令tspan=[t0,t1,t2,…,tf]（要求是单调的）。 Y0 包含初始条件的向量。 Options 用命令odeset 设置的可选积分参数。 P1,p2,… 传递给函数odefun 的可选参数。 格式 [T,Y] = solver(odefun,tspan,y0) %在区间tspan=[t0,tf]上，从t0到tf ，用初始条 件y0求解显式微分方程y’=f(t,y)。对于标量t 与列向量y ，函数f=odefun(t,y)必须返回一f(t,y)的列向量f 。解矩阵Y 中的每一行对应于返回的时间列向量T 中的一个时间点。要获得问题在其他指定时间点t0,t1,t2,…上的解，则令tspan=[t0,t1,t2,…,tf]（要求是单调的）。 [T,Y] = solver(odefun,tspan,y0,options) %用参数options （用命令odeset 生成） 设置的属性（代替了缺省的积分参数），再进行操作。常用的属性包括相对误差值RelTol （缺省值为1e-3）与绝对误差向量AbsTol （缺省值为每一元素为1e-6）。 [T,Y] =solver(odefun,tspan,y0,options,p1,p2…) 将参数p1,p2,p3,..等传递给函数 odefun ，再进行计算。若没有参数设置，则令options=[]。 1．求解具体ODE 的基本过程： （1）根据问题所属学科中的规律、定律、公式，用微分方程与初始条件进行描述。 F(y,y’,y’’,…,y (n),t) = 0 y(0)=y 0,y’(0)=y 1,…,y (n-1)(0)=y n-1 而y=[y;y(1);y(2);…,y(m-1)]，n 与m 可以不等 （2）运用数学中的变量替换：y n =y (n-1),y n-1=y (n-2),…,y 2=y 1=y ，把高阶（大于2阶）的方 程（组）写成一阶微分方程组：????????? ???=?????? ????????′′′=′)y ,t (f )y ,t (f )y ,t (f y y y y n 21n 21M M ，????????????=????????????=n 10n 210y y y )0(y )0(y )0(y y M M （3）根据（1）与（2）的结果，编写能计算导数的M-函数文件odefile 。 （4）将文件odefile 与初始条件传递给求解器Solver 中的一个，运行后就可得到ODE 的、在指定时间区间上的解列向量y （其中包含y 及不同阶的导数）。 2．求解器Solver 与方程组的关系表见表2-3。 表2-3 函数指令 含 义 函 数 含 义 ode23 普通2-3阶法解ODE odefile 包含ODE 的文件 ode23s 低阶法解刚性ODE odeset 创建、更改Solver 选项 ode23t 解适度刚性ODE 选项 odeget 读取Solver 的设置值 ode23tb 低阶法解刚性ODE odeplot ODE 的时间序列图 求解器 Solver ode45 普通4-5阶法解ODE 输出 odephas2 ODE 的二维相平面图 常微分方程初值问题数值求解 学生：张玉娟 学号：3070942232 指导老师：梁鹏 摘 要：数值求解是指在没有办法知道未知函数的解析表达式的情况下，近似计算未知函数在其定义域 中的某些离散点上的函数值。本文利用欧拉方法、梯形方法和龙格_库塔方法三种单步法基于Matlab 求解常微分方程初值问题。 关键词：常微分初值问题；数值求解；欧拉方法；梯形方法；龙格_库塔方法；Matlab 实现 1 引言 很多科学技术和工程问题常用微分方程的形式建立数学模型，因此微分方程的求解是很有意义的。建立微分方程只是解决问题的第一步，通常需要求出方程的解来说明实际现象，并加以检验。如果能得到解析形式的解固然是便于分析和应用的，但是我们知道，虽然求解常微分方程有各种各样的解析方法，但解析方法只能用来求解一些典型的方程，而对于绝大多数的微分方程问题，很难或者根本不可能得到它的解析解，实际问题终归结出来的微分方程主要靠数值解法。因此，研究微分方程求解的数值方法是非常有意义的。本文介绍了欧拉法、梯形法和四阶龙格_库塔方法三种单步法，通过Matlab 的平台运行。 2 建模 常微分方程初值问题的数学模型是：求y ()y x =，使之满足 0'()(,)()y x f x y x b y a y =≤≤??=?，(a ),, (1) 其中0y ，a ，b 是已知的常数，(,)f x y 是已知函数，且满足条件： (,)(,')'f x y f x y L y y -≤-， 式中的L 是不依赖于y ，'y 的常数，称为利普希茨（Lipschitz ）常数。 由常微分方程理论知识可知，上述问题存在唯一解()y x 。现在的目标就是计算区间[],a b 上等节点i x a ih =+处该未知函数的函数值()i y x ，其中()/h b a N =-，0,1,2,...,i N =。为此，可以用求解常微分方程问题的单步法，即欧拉方法、梯形方法和龙格_库塔方法求解。 3 算法求解和程序设计 3.1 欧拉方法 将微分方程离散化，用向前商 1()() n n y x y x h +-代替微分()n y x ，代入（1）中的微分方程，可得 1()() (,())1,2,3,...n n n n y x y x f x y x n h +-==，（） 化简可得 ┊┊┊┊┊┊┊┊┊┊┊┊┊装┊┊┊┊┊订┊┊┊┊┊线┊┊┊┊┊┊┊┊┊┊┊┊┊ 常微分方程的几种数值解法 数学与应用数学肖振华指导教师张秀艳 【摘要】自然界与工程技术中的很多现象，可以归结为微分方程定解问题。其中，常微分方程求解是微分方程的重要基础内容。但是，对于许多的微分方程，往往很难得到甚至不存在精确的解析表达式，这时候，数值解提供了一个很好的解决思路。，针对于此，本文对常微分方程数值解法进行了简单研究，主要讨论了一些常用的数值解法，如欧拉法、改进的欧拉法、Runge—Kutta方法、Adams预估校正法以及勒让德谱方法等，通过具体的算例，结合MA TLAB求解画图，初步给出了一般常微分方程数值解法的求解过程。同时，通过对各种方法的误差分析，让大家对各种方法的特点和适用范围有一个直观的感受。 【关键词】常微分方程数值解法MA TLAB 误差分析 【Abstract】 Many phenomena in nature and engineering can be attributed to the definite solution of the problem for differential equations. Among them, the ordinary differential equation solving is an important foundation for the content of the differential equations. However, many of the differential equations are often difficult to obtain accurate analytical expression .At this time, the numerical solution provides a good idea. For the Numerical Solution of Ordinary Differential Equations in this article, we focuses on some commonly used numerical solution, such as the Euler method, improved Euler method, Runge-Kutta method, Adams predictor corrector method as well as newer spectral methods. Through specific examples, combined with MATLAB solving and drawing, we initially know the solution process of general numerical solution of ordinary differential equations . At the same time, according to the error analysis of various methods , everyone has an intuitive feel of the characteristics and scope of the various methods. 【Keywords】Ordinary Differential Equations Numerical Solution MATLAB error analysis 第八章 常微分方程数值解 姓名 学号 班级 习题主要考察点：欧拉方法的构造，单步法的收敛性和稳定性的讨论，线性多步法中亚当姆斯方法的构造和讨论。 1 用改进的欧拉公式，求以下微分方程 ]1,0[1)0(2∈?????=-='x y y x y y 的数值解（取步长2.0=h ），并与精确解作比较。（改进的尤拉公式的应用） 解：原方程可转化为 x y y y 22 -='，令22 y z =，有x z dx dz 22-=- 解此一阶线性微分方程，可得 12+= x y 。 利用以下公式 )4,3,2,1,0()(21)2(2.0)2(2.01=???? ?????+=-?+=-?+=+i y y y y x y y y y x y y y c p i p i p i c i i i i p 求在节点)5,4,3,2,1(2.0=?=i i x i 处的数值解i y ，其中，初值为1,000==y x 。 MATLAB 程序如下： x(1)=0;%初值节点 y(1)=1;%初值 fprintf('x(%d)=%f,y(%d)=%f,yy(%d)=%f\n',1,x(1),1,y(1),1,y(1)); for i=1:5 yp=y(i)+0.2*(y(i)-2*x(i)/y(i));%预报值 yc=y(i)+0.2*(yp-2*x(i)/yp);%校正值 y(i+1)=(yp+yc)/2;%改进值 x(i+1)=x(i)+0.2;%节点值 yy(i+1)=sqrt(2*x(i+1)+1);%精确解 fprintf('x(%d)=%f,y(%d)=%f,yy(%d)=%f\n',i+1,x(i+1),i+1,y(i+1),i+1,yy(i+1)); end 程序运行的结果如下： x(1)=0.000000, y(1)=1.000000, yy(1)=1.000000 x(2)=0.200000, y(2)=1.220000, yy(2)=1.183216 x(3)=0.400000, y(3)=1.420452, yy(3)=1.341641 x(4)=0.600000, y(4)=1.615113, yy(4)=1.483240 x(5)=0.800000, y(5)=1.814224, yy(5)=1.612452 x(6)=1.000000, y(6)=2.027550, yy(6)=1.732051常微分方程数值解法
常微分方程数值解
常微分方程初值问题的数值解法
常微分方程的数值解法
常微分方程数值解法
MATLAB求解常微分方程数值解
常微分方程的数值解
常微分方程数值解法
常微分方程初值问题数值解法
常微分方程数值解法
(整理)常微分方程数值解法
常微分方程常用数值解法.
常微分方程数值解法
Ch5常微分方程数值解法
数值分析--第9章 常微分方程数值解
常微分方程数值解
常微分方程初值问题数值解法
常微分方程的几种数值解法
常微分方程数值解