04实验4微分方程数值解1(4)

0.8000 1.6498 1.6125
return o
0.9000 1.7178 1.6733 x1 x2 x3 x4 1.0000 1.78x48 1.7321
1.4 改进的欧拉公式
例 用改进的欧拉公式解方程
先由向前欧拉公式算出yn+1的预测值 y yny1，2yx再,把y(它0)代 1替梯形
y(1)=y0;
yn1
yn
h[ 2
f
(xn ,
yn )
f
( xn1 ,
re0turn 1.0000 yn001)..]12,00n00000,1, 121,..01985419
for n=1:m
称为改进的欧拉公式，它可写作
k1=feval(f,x(n),y(n));
0.3000 1.2662 0.4000 1.3434
yk 1 yk 1 , k 2h
1,2,, n 1
f
( x0 )
3y0
4 y1 2h
y2
,
f
( xn )
yn2
4 yn1 2h
3yn
后两个公式是由二次插值得到，目的是保持在端点处的精度O(h2) 高阶导数的近似公式一般要用插值多项式得到，下面是二阶
导f数(x的) 近((x似fx0(公axx11式)))((：xx0f(xxa22)) yh0)((x2x1
0.1000 3.9854 0.2924 15.0000 2.0000
0.1500 5.9445 0.6906 15.0000 3.0000
0.2000 7.8515 1.2899 15.0000 4.0000
0.2500 9.6705 2.1178 15.0000 5.0000
1.0000 1.7379
改进的欧拉公式可以从几何上加以解释如下：
yn1
k1
yn h k1 f (xn , yn )
2
k2
k2 f (xn1, yn hk1)
以Pn,Q两点处的斜率 之平均值为斜率，求 得Pn+1
k1 k2
k1
2
k2
Q(xn+1,yn+hk1)
y
y = f (x)
Pn(xn,yn)
明确的几何C意义：
后差公式误差为O(h)
o
a-h a a+h
x 中点公式误差O(h2)
将区间[a,
b]
n等分，步长
h
b
n
a
.
当函数 y= f (x)在分点
上用离散数值表示为( xk , yk)，a=x0< x1 <…< xn = b 时, 函数
在分点的导数值由中点公式得到三点公式：
f
( xk )
end
y
0.3000 1.2774 1.2649
z=[x' y‘];
p2
p3
0.4000 1.3582 1.3416 p40.5000 1.4351 1.4142
0.6000 1.5090 1.4832
function dy=odpe1121(x,y)
0.7000 1.5803 1.5492
dy=y-2*x/y;
人口(百万) 76.0 92.0 106.5 123.2 131.7 150.7 179.3 204.0 226.5 251.4 281.4
若记时刻 t 的人口为x(t), 则人口的相对增长率 r(t)=x’(t)/x(t) 则由三点公式计算得年增长率为
表2 美国人口的年增长率 年 份 1900 1910 1920 1930 1940 1950 1960 1970 1980 1990 2000 年增长率r(%) 2.2 1.66 1.46 1.02 1.04 1.58 1.49 1.16 1.05 1.04 1.16
es=nky[d2(xn'=+yfe1'];v)=akyykl(n12(fn,1x)(+ffnh(y(+*xxn1(nnk),1,1yy,h2+ny()(nknk)21+)h/h2kk*;12k))1);
0.5000 0.6000 0.7000 0.8000 0.9000
1.4164 1.4860 1.5525 1.6165 1.6782
导数有以下三个近似公式：
f (a) f (a h) f (a) h
f (a) f (a) f (a h) h
y
f (a) f (Ta h) f (a h) 以上三个公式分别
A
B
y=2hf(x)
表示了线段AB、CA
它们依次称为前差公式、后差公式和中点公式。这些公式都
和CB的斜率。前、
4.3.3 龙格-库塔方法
注意到 y(xn1) y(xn ) hy(xn h)
选取 xn , xn a h 两个点上的导数作线性组合
yn1 k1
yn h(l1k1
f (xn , yn )
l2k2 )
k2 f (xn ah, yn ahk1)
这里l1,l2为加权系数
a (0< a<1),3个数均为
Pn+1(xn+1,yn+1)
o
xn
xn+1
x
欧常拉用方的法方可法以是推:广令到y解1 高y阶, y方2 程 和y,方y3程组y，,如,：yn y(n1).
zynn向11 转前yyyy13z2n换的ynn为欧yyy243ahdd一d拉dh1yxgxzf(阶公x(()xxg常y式fnn(n(,,xx微为yy,, nyyan分：,,,,2zzz(z))线xnn,)))y性y(nnx方01)程0,y1组0,, 2z:a(,nx0方1)(zyx程nn)z11y0组2 的a向nzyn(nx量)y形1h式gff
f (xn , yn )
l2k2 )
k2 f (xn ah, yn ahk1)
称为二阶龙格-库塔公式。
（39）
可以证明，当在小区间[ xn, xn+1]内只取2 点的二 阶龙格-库 塔公式的最高精度为2 阶（O(h3)）。
类似地，在小区间选四 个点，适当选取加权系数， 得到精度为四阶O(h5)的四阶 龙格-库塔公式：
公式的右端的yn+1，作为校正，即
yf=uandcvteiounla(d'oyd=eo1d2e11'2,01,(1x,,1y,)10)
function s=adveula(f,a,b,y0,h)
yd=y=y-2*x/y;
x=a:h:b;my=nf1looyr(n(bh-af ()/xhn ,);yn )
速a度=2前0;往b=拦40截;c=。15用; 雷达进行跟踪时，可保持缉私艇的速度方向始 终模sd指型=xs=向，q[r(tc走确(-(xc私定(-1x船缉)()1/)。私s);^(建艇a2*+t立追(-ax*任上(2t-)意走x)(/s2时私])*)刻船^b2; 缉的)%; 私位以艇置列的，向位求量置出形和追式缉上表私的示艇时方航间程线。的数学
自然界中不同物种之间存在着一种非常有趣的相互依存相互制约的生存方式种群甲靠丰富的自然资源生长而种群乙靠捕食甲为生鱼和鲨鱼兔和山猫落叶松和蚜虫都是这种生存方式的典型
实验4
4.2 数值微分
一、数值微分公式 当函数以离散数值形式给出时，用离散方法 近似地计算函数y= f (x)在某点的导数值，这就是数值微分。
模再型编设写在一t=个0时脚刻本缉M私文艇发件现走 Y
私c船lea，r 走私船以速度a平行于Y轴
正ts向=0行:0驶.05，:0缉.5私; 艇按速度b以指
向x走0=私[0船0]的'; 方向行驶(b>a)，在任
Q(c, at)
意[时t x刻]=otd缉e45私(@艇ji位si,于ts,点x0P);(x,y),则 M文件以ex4t1作为文件
4.3 欧拉方法和龙格-库塔(Runge-Kutta)方法
设有常微分方程初值问题：
y' f (x, y)
常<1欧拉方法常在x微唯微2 <分一向向梯改分…方解方前后形进程<.y程已 试 欧欧公欧nx数结n知 1拉 拉 式 拉求h<值论…y::解:y上(fn：x(y适yyy的x0yy(n计nn0nn当)近x1,111算1n光y似)0（f滑)yy值yy,ynn(精n0，yynxy，确'n,(h(对h2(1通xfxy）[f0y(0)ff(f常x(解满()(nx(xxfx,选xnn)y足n,y,(,)fn(yhn取xy1yL),x()))nnix在,y相n)p)n,nsx一等1cfy)hf(,h)2系的ix(f1ty[,xz,(x列y计n条n2xn21,)离,n算件,1.xy.)1.散n步,n,h则L1y1点长)(]方]nyn11xh程)0。<存yh02,x1)1, 2,
同样可推广到解微分方 程组和ቤተ መጻሕፍቲ ባይዱ阶方程。
yn1
yn
h 6
(k1
2k2
2k3
k4 )
k1
k2
k3
f (xn , yn )
f
(
xn
h 2
,
yn
f
( xn
h 2
,
yn
hk1 ) 2 hk2 ) 2
k4 f ( xn h, yn hk3 )
4.4 龙格-库塔方法的MATLAB实现
先对所解方程或方程组（高阶方程）编写ODE文件（M 函数文件，Ordinary Differential Equation）保存在MATLAB 设置的搜索路径，然后在解方程的脚本文件中调用：
4.3.1 欧拉方法 1例.1 用向向前前欧欧拉拉公公式式解方程
y> >zy=fyo2nyexu1,hlya((0y'o)n de11f2(1x'n,0, ,y1n,)1,0.1);
>>y1=sqrt(1+2*x);
funct设io初n 值z=为foeyu(xl0a)(=f,ay,0b，,y0yn,+h1) = y>n>+sohlf=([zxny,1y’n] ), n = 0, 1, 2, …
待定常数。
l1,l2 ，a 的选取标准：让精度更高
由二元函数的泰勒展式（推导在教材P74）,得
当满足l1 l2=1, 2a l2=1 时，精度达到2阶。显然解不唯一。 特别当l1 l21/ 2，a 1时就是改进的欧拉公式。
满足l1 l2=1, 2a l2=1 的公式
yn1 k1
yn h(l1k1
O
cX
gtdextt('x');g(tecxt('yx') 2 (at y)2
设走私船在t时刻的位置为x1(t),y1(t),则
x1(t) c 15, y1(t) at 20t
时刻t 缉私艇的位置 走私船的位置
0
0
0 15.0000
0
0.0500 1.9984 0.0698 15.0000 1.0000
[t, x]=ode23(@f,ts,x0,options) 或
[t, x]=ode45('f',tspan,x0,options)
ts=[a, b]表示自变量的 初值到终值的区间
或ts=[t0,t1,t2,…tf]
options用于设置误差限（缺省时相对误差为10-3，绝对误差 为10-6），用语句实现为：
[tdxx] b cosa , dy b sin a plodtt(t,x);grid dt
即 gfitgddeuxxtrte(';x(t)')(;cgtexxtb()'(y2c(t)')(;xpa)atusey)2
pldoty(x(:,1),x(:,2)b);(garitd y)
名径存中P盘。(x,，y)保存a在搜索R路(c, y)
(xn , (xn , (x)
yn yn
, ,
zn zn
) )
其中x0 , y0 , z0已知
解线性方程组,可采用向量函数表示.
对于高阶方程,需先降阶化为一阶常微分方程组,然后再按上方
法求解.如对高阶线性方程
y(n) a1(x) y(n1) a2 (x) y(n2) an1(x) y ' an (x) y f (x)
fh(2xxa00)))((xx1f x(x2a2))y1h)
(
(x x2
x0 x0
)( )(
x x1) x2 x1)
y2
diff(x): [x2-x1 , x3-x2,… xn-xn-1]
二、应用实例——人口增长率 已知20世纪美国人口统计数据如下，试计算表1 中这些年
份的人口增长率。
表1 20世纪美国人口统计数据 年份 1900 1910 1920 1930 1940 1950 1960 1970 1980 1990 2000
x=a:h:b;m=floor((b-a)/h);
sol=
y这(1就)=叫y0向; 前欧拉公式。其精度为1 阶0 （局部截1.断00误00差为1O.00(h020))。
for n=1:m
0.1000 1.1000 1.0954
y(n+1)=y(n)+h*feval(f,x(n),y(n)); 0.2000 1.1918 1.1832
options=odeset('reltol',rt,'abstol',at)
§设 4a=.250,b=应40,用c=1问5; 题求解
例有先f1一u、n建艘c海t立i走o上n一私缉d船x个私=正j函is以海i(数t,一防xM)定某文速部件度缉向私正艇北上方的向雷行达驶发，现缉正私东艇方立向即C海以里最处大