第十章 常微分方程(组)求解(7)

第三篇 第十章 常微分方程（组）求解

上的数值解','y/dx=f(x,y)，y(x0)=y0的精确解y=f(x)')

end

%计算误差. for k=2:n+1

wucha(k)=norm(Y(k-1)-Y(k)); wch(k)=norm(fxy(k)-Y(k)); end

X=X(1:k); Y=Y(1:k,:);fxy=fxy(1:k,:); n=1:k;wucha=wucha(1:k,:); wch=wch(1:k,:);

P=[n',X,Y,fxy,wch,wucha];

用常用四阶龙格-库塔方法求解常微分方程初值问题（10.5）的数值解的MATLAB主程序2

function [x,y]=RKc4(x0,b,y0,h) for x=x0:h:b

k1=funfcn(x0,y0);

k2=funfcn(x0+h/2,y0+h*k1/2); k3=funfcn(x0+h/2,y0+h*k2/2);

k4=funfcn(x0+h,y0+h*k3);

y=y0+h*(k1+2*k2+2*k3+k4)/6;

x=x0+h;[x,y], x0=x;y0=y; end [x,y];

例10.5.3 用常用的四阶龙格-库塔公式求初值问题

2xy?dy,0?x?2,??1?2 dx1?x???yx?0?0的数值解，h=1/4,并计算与精确解的误差，在同一图形窗口画出例10.5.1、例10.5.2和

例10.5.3精确解和数值解的图形.

解 方法1

(1) 编写并保存名为funfcn.m和fun.m的MATLAB程序 (2)在MATLAB工作窗口输入下面的程序

>> x0=0;b=2; y0=0;h=1/4; subplot(3,1,1)

C=[1/4,3/4,2/3,2/3];

[k,X,Y,fxy,wch,wucha,P]=RK2(@funfcn,@fun,x0,b,C,y0,h) subplot(3,1,2)

c1=1/6;c2=4/6;c3=1/6;a2=1/2; a3=1;b21=1/2;b31=-1;b32=2;

C=[c1,c2,c3,a2, a3,b21,b31,b32];

[k,X,Y,fxy,wch,wucha,P]=RK3(@funfcn,@fun,x0,b,C,y0,h) subplot(3,1,3)

c1=1/6;c2=2/6;c3=2/6;c4=1/6;a2=1/2; a3=1/2; a4=1;b21=1/2;b31=0;b32=1/2; b41=0;b42=0;b43=1;

C=[c1,c2,c3, c4,a2, a3, a4,b21,b31,b32,b41,b42,b43]; [k,X,Y,fxy,wch,wucha,P]=RK4(@funfcn,@fun,x0,b,C,y0,h)

将运行后计算的结果（略），画出精确解和数值解的图形.

方法2 在MATLAB工作窗口输入下面的程序

>> x0=0;b=2; y0=0;h=1/4; [x,y]=RKc4(x0,b,y0,h)

运行后输出结果（略）.

10.5.5 自适应龙格-库塔方法及其MATLAB程序

200．

第三篇 第十章 常微分方程（组）求解

(一) 常微分方程（组）的刚性和非刚性 我们先对刚性方程给直观的描述。

自然界的各种现象发生的每一过程都是及其复杂的，往往包含许多子过程及其它们之间的相互作用，其中有些子过程表现为快的变化，而有些则相对慢些，且变化的速度可以相差非常大的数量级，相应地描述这些子过程的常微分方程（组）的解中也将包含快变化和慢变化的分量，如果在一个过程中的快变量过程与慢变量过程的变化速度相差非常大，在数学上称这种过程具有“刚性”，而描述这种过程的常微分方程（组）称为刚性方程（组），也称为病态方程（组）或坏条件方程（组）。否则称为非刚性方程（组）。

??Ay?f，对于一阶常系数线性微分方程组y如果它的雅可比矩阵An?n的特征值相差十

分悬殊，则称此微分方程组是刚性的。“刚性”问题在控制系统工程、电子网络、生物学，物理学及化学动力学过程中经常遇到

考虑线性微分方程

du?Au?f(x), dx其中A是n阶常数矩阵，u和f是n维向量函数，

?u1(x)??f1(x)??u(x)??f(x)?2?,f(x)??2?. (I) u????????????u(x)f(x)?n??n?记矩阵A的特征值为?i(i?1,2,...,n)，相应的特征向量为zi，则上述方程组的通解是 u??ce?ii?1nixzi??(x),

其中?(x)是方程组的特解，常数c1,c2,...,cn可以由初始条件确定。如果

Re(?i)?0(i?1,2,...,n)，则当x???时，?ciezi?0。称?cie?ixzi是方程组

?ixi?1i?1nn的渐态解，?(x)是稳态解。用数值方法求方程组的稳态解上四需要计算到渐态解中变化最慢的分量可以忽略为止。

定义1：对于常系数线性微分方程组(I)，若矩阵A的特征值?i(i?1,2,...,n)满足 （1）Re(?i)?0;

201．

第三篇 第十章 常微分方程（组）求解

（2）s?maxRe(?i)iminRe(?i)i??1,

则称方程组(I)是刚性(stiff)方程组，s称为刚性比。

通常当刚性比s达到O(10p)(p?1)时，就认为微分方程组是刚性的。s的量级为10，称作临界刚性，s越大刚性越严重。

?dui?fi(x,u1,u2,...,un),x?x?b,?对于一般的方程组?dx (II)

??ui(x0)?ui0,i?1,2,...,n,其中ui是未知函数，ui0是初值。其Jacobi（雅可比）矩阵

??f1??u?1??f2?uJ(x)???1?...???fn???u1?f1?u2?f2?u2...?fn?u2?f1??un???f2?...?un?

?......???fn?...?un??...其中偏导数在在方程组(II)的真解处取值，它们都是x的函数。记

?i(x)(i?1,2,...,n)是矩阵J(x)的特征值。

定义2：对于微分方程组(II)，若对x?[x0,b]成立： （1）Re(?i(x))?0(i?1,2,...,n);

（2）s?maxRe(?i)iminRe(?i)i??1,

则称方程组(II)是刚性(stiff)方程组，s称为在刚性比。

刚性方程组又称为病态(ill?conditioned)方程组，在化学反应、自动控制、电力系统等方面常遇到这种问题。对于刚性问题应当采用绝对稳定性比较好的数值方法求解，

202．

第三篇 第十章 常微分方程（组）求解

以至于步长可以取的大一些，最好是对步长不加限制。 (二) 解常微分方程初值问题的MATLAB库函数

表10-12 解常微分方程初值问题的库函数

odeij 问题类型 精度 适用对象

ode45 非刚性 中等 多数情况下可优先适用，但不能用来解刚性问题 ode113 非刚性 低到高 不能解刚性问题，当误差容限要求严格时，结果

较ode45 好。 ode23 非刚性 较低 可用来解中等刚性问题，或误差允许范围较宽的

问题。 ode15s 刚性 低到高 可用来解刚性问题,当采用ode45失败或效果很

差时，可考虑适用。 ode23s 刚性 低 可用来解刚性问题，当误差容限要求较宽时，结

果较ode15s好。 ode23tb 刚性 低 可用来解刚性问题，当误差容限要求较宽时，结

果较ode15s好。 ode23t 刚性 低 可用来解刚性问题，但要求无数值衰减。

调用格式一：[T,Y] = ode ij (@ODEFUN,TSPAN,Y0)

调用格式二：[T,Y] = ode ij (@ODEFUN,TSPAN,Y0,OPTIONS)

调用格式三：[T,Y] = ode ij (@ODEFUN,TSPAN,Y0,OPTIONS,P1,P2,...) 调用格式四:[T,Y,TE,YE,IE]=odeij(@ODEFUN,TSPAN,Y0,OPTIONS,P1,P2,...)

?dy2x?,0?x?1.2,?例10.5.4 分别用二阶数值方法求初值问题?dx3y2的数值解，

?yx?0?1??4精确到10,计算它与精确值的绝对误差和相对误差，并画出精确解和数值解的图形.

解 （1）先求精确解，输入程序

>> y=dsolve('(Dy)-2*x/(3*y^2)=0','y(0)=1','x') (2) 编写并保存名为funfcn.m的MATLAB程序如下

function f=funfcn(x,y)

f=2*x/(3*y^2);

(3)在MATLAB工作窗口输入下面的程序

>> options=odeset('RelTol', 1e-4,'AbsTol',1e-4); [t,y]=ode23(@funfcn,[0 1.2],1,options) , yf=(t.^2+1).^(1/3)

plot(t,y(:,1),'rp',t,yf,'bo'); grid

xlabel('自变量 X'), ylabel('因变量 Y')

legend('用二阶龙格-库塔方法计算dy/dx=2x/(3y^2)，y(0)=1在

[0,1.2]上的数值解',' dy/dx=2x/(3y^2)，y(0)=1的精确解y=f(x)')

juew=yf(:,1)-y(:,1), xiangw=juew./yf(:,1),

203．

第三篇 第十章 常微分方程（组）求解

[t,y,yf,juew, xiangw] 运行后计算的结果（略），并画出精确解和数值解的图形.

例10.5.5 对例10.5.4也可以将ode23函数换成ode45函数求其初值问题的数值解，精确到10时得到的结果是相同的。试对这两种方法进行比较。 解：编写名为tode23.m的M文件 function tode23 tic;p1=flops;

options=odeset('RelTol',1e-4,'AbsTol',1e-4); [x,y]=ode23(@function,[0 1.2],1,options); p2=flops;tode23=toc,pode23=p2-p1 在MATLAB工作窗口输入文件名 >> tode23

运行后输出运行时间如下 tode23=

0.2800

把名为tode23.m的M文件中的ode23改成ode45，运行文件名 >> tode45 后，得

tode45= 0.2200

从此例所得结果可以看出，ode23比ode45确实慢些，但是ode23比方法要快。

例10.5.6 先判别方程组

?4?dy?dx??2y?z?2sinx,??dz?998y?999z?999(cosx?sinx), ?dx?yx?0?2??zx?0?3?是否是刚性的.再分别用ode45函数和ode15s函数求在[0, 1.2]上的数值解，精确到10,

将计算结果与精确值比较，并画出精确解和数值解的图形.

解 （1）在MATLAB中求刚性比，输入程序

>> A=[-2 1;998 -999];Tez1=max(abs(real(eig(A)))); Tez2=min(abs(real(eig(A))));Gxb=Tez1/Tez2 运行后输出刚性比为 Gxb=

1000

因为刚性1000比远大于1，所以原方程组为刚性方程组。 (2)求精确解。输入程序 >> syms x y z

f='(Dy)-2*sin(x)-z+2*y=0,(Dz)+999*sin(x)-999*cos(x)+999

*z-998*y=0';

[y,z]=dsolve(f,'y(0)=2,z(0)=3','x')

运行后屏幕显示常微分方程组在给定初始条件下的精确解y和z如下

y= z=

2*exp(-x)+sin(x) 2*exp(-x)+cos(x) (3)分别用函数ode45和ode15s函数解刚性方程的数值解，误差10画出精确解和数值解的图形。

(1)编写并保存名为funfcn.m的MATLAB程序如下 function f=funfcn(x,u)

f=[-2 1;998 -999]*u+[2*sin(x);-999sin(x)+999*cos(x)];

(2)因为原方程组为刚性方程组，所以在MATLAB工作窗口输入下面的程序 function li1056

204．

?1?1

百度搜索“77cn”或“免费范文网”即可找到本站免费阅读全部范文。收藏本站方便下次阅读，免费范文网，提供经典小说综合文库第十章 常微分方程(组)求解(7)在线全文阅读。