ODE45竟然不如改进欧拉法及龙格库塔法—有程序验证

simber

我做了个数值试验，研究各种数值解法的长期行为，分别采用了四种途径计算椭圆的微分方程组：
1：欧拉算法  2：改进欧拉法   3：四阶龙格库塔法   4：matlab的ODE45
试验步长都取一样，绕椭圆画了几圈，发现而改进的欧拉算法和四阶标准龙格－库塔法稳定性都很好，基本没有偏离原椭圆轨道，同时意外发现matlab自带的ODE45比改进的欧拉算法和四阶标准龙格－库塔法稳定性都要差，从局部放大图可以明显看出。。。


我现在都不相信ODE45了，用自己编的龙格－库塔算。。。请问高手是不是ODE45有什么设置不对？？能通过设置改进精度吗？

附程序，共5个m文件。
★文件“jsff04.m”
clear all ;clc;close all
n=1e4;
h=2e-3;
a=1;b=2;
yy1=f1(a,b,h,n);
yy2=f2(a,b,h,n);
yy3=ronge(a,b,h,n);
figure('Name','ff1','Position',[1 1 680 400])  %作图设置
subplot(1,2,1)
plot(yy1(1,:),yy1(2,:))
axis([-2,2,-2.5,2.5]);axis equal
title('欧拉算法')
subplot(1,2,2)
plot(yy1(1,:),yy1(2,:))
axis([0.2,0.3,1.96,1.968])
title('欧拉法局部放大')
saveas(gcf, '欧拉算法.emf')  
figure('Name','ff2','Position',[1 1 680 400])  %作图设置
subplot(1,2,1)
plot(yy2(1,:),yy2(2,:))
axis([-2,2,-2.5,2.5]);axis equal
title('改进欧拉算法')
subplot(1,2,2)
plot(yy2(1,:),yy2(2,:))
axis([0.2,0.3,1.96,1.968])
title('改进欧拉算法局部放大')
saveas(gcf, '改进欧拉算法.emf')  
figure('Name','ff3','Position',[1 1 680 400])  %作图设置
subplot(1,2,1)
plot(yy3(1,:),yy3(2,:))
axis([-2,2,-2.5,2.5]);axis equal
title('龙格库塔算法')
subplot(1,2,2)
plot(yy3(1,:),yy3(2,:))
axis([0.2,0.3,1.96,1.968])
title('龙格库塔算法局部放大')
saveas(gcf, '龙格库塔算法.emf')  
figure('Name','ff4','Position',[1 1 680 400])  %作图设置
subplot(1,2,1)
[t,y]=ode45('fode',[h:h:h*n],[0,2],[],a,b);
plot(y(:,1),y(:,2))
axis([-2,2,-2.5,2.5]);axis equal
title('Matlab ode45')
subplot(1,2,2)
plot(y(:,1),y(:,2))
axis([0.2,0.3,1.96,1.968])
title('Matlab ode45结果局部放大')
saveas(gcf, 'Matlab ode45.emf')  
★文件“f1.m”
function y=f1(a,b,h,n)
x(1)=0;
y(1)=b;
for i=1:n
    y(i+1)=y(i)-h*b*x(i);
    x(i+1)=x(i)+h*a*y(i);
end
y=[x;y];
★ 文件“f2.m”
function y=f2(a,b,h,n)
x(1)=0;
xx(1)=0;
y(1)=b;
yy(1)=b;
for i=1:n
    xx(i+1)=x(i)+h*a*y(i);
    yy(i+1)=y(i)+h*(-b*x(i));
    x(i+1)=x(i)+h/2*(a*y(i)+a*yy(i+1));
   y(i+1)=y(i)+h/2*((-b*x(i))+(-b*xx(i+1)));
end
y=[x;y];
★文件“ronge.m”
function y=ronge(a,b,h,n)
x(1)=0;
y(1)=b;
for i=1:n
    Kx1=a*y(i);
    Ky1=-b*x(i);
    xx1=x(i)+(h/2)*Kx1;
    yy1=y(i)+(h/2)*Ky1;
    Kx2=a*yy1;
    Ky2=-b*xx1;
    xx2=x(i)+(h/2)*Kx2;
    yy2=y(i)+(h/2)*Ky2;
    Kx3=a*yy2;
    Ky3=-b*xx2;
    xx3=x(i)+(h)*Kx3;
    yy3=y(i)+(h)*Ky3;
    Kx4=a*yy3;
    Ky4=-b*xx3;
x(i+1)=x(i)+(h/6)*(Kx1+2*Kx2+2*Kx3+Kx4);
    y(i+1)=y(i)+(h/6)*(Ky1+2*Ky2+2*Ky3+Ky4);
end
y=[x;y];
★文件“fode.m”
function ydot=fode(t,y,flag,a,b)
ydot=[1*y(2);-2*y(1)];

sogooda

这个问题用odeset改下精度应该就解决了。

ChaChing

个人欣赏楼主这种态度! 可惜个人时间/能力有限

sogooda

clear;clc
n=1e4;
h=2e-3;
a=1;b=2;
yy1=f1(a,b,h,n);
yy2=f2(a,b,h,n);
yy3=ronge(a,b,h,n);
figure('Name','ff1','Position',[1 1 680 400])  %作图设置
subplot(1,2,1)
plot(yy1(1,:),yy1(2,:))
axis([-2,2,-2.5,2.5]);axis equal
title('欧拉算法')
subplot(1,2,2)
plot(yy1(1,:),yy1(2,:))
axis([0.2,0.3,1.96,1.968])
title('欧拉法局部放大')
saveas(gcf, '欧拉算法.emf')  
figure('Name','ff2','Position',[1 1 680 400])  %作图设置
subplot(1,2,1)
plot(yy2(1,:),yy2(2,:))
axis([-2,2,-2.5,2.5]);axis equal
title('改进欧拉算法')
subplot(1,2,2)
plot(yy2(1,:),yy2(2,:))
axis([0.2,0.3,1.96,1.968])
title('改进欧拉算法局部放大')
saveas(gcf, '改进欧拉算法.emf')  
figure('Name','ff3','Position',[1 1 680 400])  %作图设置
subplot(1,2,1)
plot(yy3(1,:),yy3(2,:))
axis([-2,2,-2.5,2.5]);axis equal
title('龙格库塔算法')
subplot(1,2,2)
plot(yy3(1,:),yy3(2,:))
axis([0.2,0.3,1.96,1.968])
title('龙格库塔算法局部放大')
saveas(gcf, '龙格库塔算法.emf')  
figure('Name','ff4','Position',[1 1 680 400])  %作图设置
subplot(1,2,1)

options=odeset('RelTol',1e-9);
[t,y]=ode45('fode',[h:h:h*n],[0,2],options,a,b);
plot(y(:,1),y(:,2))
axis([-2,2,-2.5,2.5]);axis equal
title('Matlab ode45')
subplot(1,2,2)
plot(y(:,1),y(:,2))
axis([0.2,0.3,1.96,1.968])
title('Matlab ode45结果局部放大')
saveas(gcf, 'Matlab ode45.emf')

simber

刚看了帮助，odeset原来还有那么多设置

sososwim

我认为ode45正确称呼应该是RKF45，即变步长的runge-kutta-fehlberg4阶和5阶算法，步长是自动调整的，不知道楼主怎么可能设置步长，该方法比前三种在工程应用上更加具有普适性，因为前三种的步长是一定的，但是并没有估计每一步的误差量，而ode45通过4阶和5阶的差值估计误差，再通过miline device调整步长，减半或者加倍（但是是有上限的）。ode45是单步法，如果要精读更高的可以用ode113，是ABM PECE线性多步法，即ADAMS-Bashful-Moulton prediction-correction （预测校正公式），但是这些算法的缺点都是只能计算连续系统的常微分方程，对于不连续系统需要修改程序，判断不连续的点。

sleepinglion

我也感觉自己编写的龙格库塔法似乎比ODE45更好！因为之前我用它们求一个系统的动力学响应时，发现自己编写的比ODE45得到的曲线更光滑！不知道为什么？我现在越来越不信赖MATLAB了，感觉它有问题！有时候一个变量符号错了，前几次运行的时候它不提示错误，然后后面运行时会突然提示错误；不一样的矩阵赋值方式，得到的矩阵结果是一样的，但计算结果却不一样！这些让我很郁闷！