matlab数值分析程序--高等数学,数值代数的matlab实现-文字版, matlab电子书, 和matlab 有关的电子书:

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


10.5.5  自适应龙格-库塔方法及其MATLAB程序
 (一) 解常微分方程初值问题的MATLAB库函数
调用格式一:[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,...)
例10.5.4  分别用二阶数值方法求初值问题的数值解,精确到,计算它与精确值的绝对误差和相对误差,并画出精确解和数值解的图形.
 
解  在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),
[t,y,yf,juew, xiangw]
运行后计算的结果(略),并画出精确解和数值解的图形.
 
例10.5.6   先判别方程组

是否是刚性的.再分别用ode45函数和ode15s函数求在[0, 1.2]上的数值解,精确到,将计算结果与精确值比较,并画出精确解和数值解的图形.
解  因为原方程组为刚性方程组,所以在MATLAB工作窗口输入下面的程序
options=odeset('RelTol', 1e-1,'AbsTol',[1e-1,1e-1]);
[t,y1]= ode15s(@funfcn,[0 10],[2 3],options)
yf=2*exp(-t)+sin(t);
zf=2*exp(-t)+cos(t);
plot(t,y1(:,1),'mo',t,y1(:,2),'rp');
hold on
plot(t,yf,'b-',t,zf,'g-.')
hold off
grid
xlabel('自变量 X'), ylabel('因变量 Y')
legend( '用ode15s函数解刚性方程的y=f(x)数值解','用ode15s函数解刚性方程的z=g(x)数值解','精确解y=f(x)','精确解z=g(x)')
运行后屏幕显示计算的结果(略),并画出精确解和用ode15s求的数值解的图形.
下面计算时间,编写并保存tode.m 的程序
function tode
tic; p1=flops;
options=odeset('RelTol', 1e-1,'AbsTol',[1e-1,1e-1]);
[t,y1]= ode45(@funfcn,[0 10],[2 3],options);
p2=flops;
tode45=toc,
pode45=p2-p1,
tic;,p3=flops;
options=odeset('RelTol', 1e-1,'AbsTol',[1e-1,1e-1]);
[t,y2]= ode15s(@funfcn,[0 10],[2 3],options);
p4=flops;
tode15s=toc,
pode15s=p4-p3
输入命令
>> tode
运行后屏幕显示计算结果
tode45 =            tode15s =
4.5100                  0.0500
由此可以看出,对于刚性的方程组,用ode15s求数值解比用ode45速度快许多.
欢迎转载,转载请注明来自一手册:http://yishouce.com/book/3/3100505.html
友情链接It题库(ittiku.com)| 版权归yishouce.com所有| 友链等可联系 admin#yishouce.com|粤ICP备16001685号-1