四阶5级龙格-库塔法出错：y从第四列值全部为nan

navycui

运行后 从第四列值全部为 nan,请高手指点

function y=lorenzeq(t,x)
sigma=10; b=8/3; r=28;
y=[-sigma*(x(1)-x(2));
     -x(1)*x(3)+r*x(1)-x(2);
     x(1)*x(2)-b*x(3)];

t_0=0 ; t_final=100; x0=[0.01;0.01;0];
t=[t_0:0.1:t_final];
n=length(t);
y=zeros(length(x0),length(t));
y(:,1)=x0;
h=0.1
for i=1:(n-1);
k1 = lorenzeq(t(i),y(:,i));
k2 = lorenzeq(t(i)+h,y(:,i)+1/4*k1);
k3 = lorenzeq(t(i)+h,y(:,i)+3/32*k1+9/32*k2);
k4 = lorenzeq(t(i)+h,y(:,i)+1932/2197*k1-7200/2197*k2+7296/2197*k3);
k5 = lorenzeq(t(i)+h,y(:,i)+439/216*k1-8*k2+3680/513*k3-845/4104*k4);
k6 = lorenzeq(t(i)+h,y(:,i)-8/27*k1+2*k2-3544/2565*k3+1859/4104*k4-11/40*k5);
y(:,i+1)=y(:,i)+h*(16/135*k1+0*k2+6656/12825*k3+28561/56430*k4-9/50*k5+2/55*k6);
end
plot(t,y)
figure;
plot3(y(:,1),y(:,2),y(:,3));
grid

[ 本帖最后由 eight 于 2007-6-8 23:10 编辑 ]

studyboy

建议用ode45求解，你自己编的程序很乱，而且执行效率很低。

navycui

这个是编程作业,不能用ode45求解

yufeng

我学的龙格法怎么和你的不一样

xjzuo

LZ的程序本身就不对, 这应该不能叫RK法．
还有，＂四阶5级龙格-库塔法＂的说法还是第一次听说．．．

bainhome

楼主关于四阶五级的说法是成立的，所编写程序是RK法的一个改进。
基本算法是分别用h和h/2两个步长两次求解，缺点是步长较小计算量大；另外，结果不好时要重新计算。
RKF45即：Runge-Kutta-Fehlbrg法，一般称为4阶五级或者四阶半RK法，目的是改善原来RK4的这一缺点，给一流程来判断是否使用了正确的步长h，每一步中用不同方法求近似解比较结果。，所以有六个数值。其实这个问题我感觉难点不在rk本身，而在于最佳步长的推导上——它是变步长的求解方式。这其实就是ode45的算法。

xjzuo

原帖由 bainhome 于 2007-6-9 12:17 发表

                               
登录/注册后可看大图

楼主关于四阶五级的说法是成立的，所编写程序是RK法的一个改进。
基本算法是分别用h和h/2两个步长两次求解，缺点是步长较小计算量大；另外，结果不好时要重新计算。
RKF45即：Runge-Kutta-Fehlbrg法，一般称为 ...

那看来是我孤陋寡闻了. 有时间我再看看这方面的程序和介绍.
------不过上面程序的计算结果的确是相当糟糕的,不知bainhome兄有无更好的建议？
(虽然用ode45可以轻易解决该问题)

[ 本帖最后由 xjzuo 于 2007-6-9 13:19 编辑 ]

bainhome

定步长RK4(自编):

1. function Y=RungeKutta4(f,xn,y0)
   
2. % xn=0:.1:1;
   
3. % y0=1;
   
4. % y_n=[];
   
5. % f=@(X1,Y1) Y1-2*X1/Y1;
   
6. y_n=[];
   
7. h=diff(xn(1:2));
   
8. for i=1:length(xn)-1
   
9.     K1=f(xn(i),y0);
   
10.     K2=f(xn(i)+h/2,y0+h*K1/2);
    
11.     K3=f(xn(i)+h/2,y0+h*K2/2);
    
12.     K4=f(xn(i)+h,y0+h*K3);
    
13.     y_n=[y_n;y0+h/6*(K1+2*K2+2*K3+K4)];
    
14.     y0=y_n(end);
    
15. end
    
16. Y=y_n;

复制代码

自适应步长RKF45(From J.K.Mathews，我做了小幅度修改，去掉了我认为不必要的测试语句，加了些注释):

1. function R=RungeKutta45(f,a,b,ya,m,tol)
   
2. % Input    - f  function
   
3. %             - a  left  endpoint of [a,b]
   
4. %             -b right endpoint of [a,b]
   
5. %             - ya initial value
   
6. %             -m initial guess for number of steps
   
7. % Output - T  solution: vector of abscissas
   
8. %             - Y  solution: vector of ordinates
   
9. 
   
10. %  NUMERICAL METHODS: Matlab Programs
    
11. % (c) 2004 by John H. Mathews and Kurtis D. Fink
    
12. %  Complementary Software to accompany the textbook:
    
13. %  NUMERICAL METHODS: Using Matlab, Fourth Edition
    
14. 
    
15. a2 = 1/4; b2 = 1/4; a3 = 3/8; b3 = 3/32; c3 = 9/32; a4 = 12/13;
    
16. b4 = 1932/2197; c4 = -7200/2197; d4 = 7296/2197; a5 = 1;
    
17. b5 = 439/216; c5 = -8; d5 = 3680/513; e5 = -845/4104; a6 = 1/2;
    
18. b6 = -8/27; c6 = 2; d6 = -3544/2565; e6 = 1859/4104; f6 = -11/40;
    
19. r1 = 1/360; r3 = -128/4275; r4 = -2197/75240; r5 = 1/50;
    
20. r6 = 2/55; n1 = 25/216; n3 = 1408/2565; n4 = 2197/4104; n5 = -1/5;
    
21. big = 1e15;
    
22. h = (b-a)/m;
    
23. hmin = h/64;% 步长自适应范围下限
    
24. hmax = 64*h;% 步长自适应范围上限
    
25. max1 = 200;% 迭代次数上限
    
26. Y(1) = ya;
    
27. T(1) = a;
    
28. j = 1;
    
29. tj = T(1);
    
30. br = b - 0.00001*abs(b);
    
31. while (T(j)<b),
    
32.   if ((T(j)+h)>br), h = b - T(j); end
    
33.   tj = T(j);
    
34.   yj = Y(j);
    
35.   y1 = yj;
    
36.   k1 = h*f(tj,y1);
    
37.   y2 = yj+b2*k1;                         if big<abs(y2) return, end
    
38.   k2 = h*f(tj+a2*h,y2);
    
39.   y3 = yj+b3*k1+c3*k2;                   if big<abs(y3) return, end
    
40.   k3 = h*f(tj+a3*h,y3);
    
41.   y4 = yj+b4*k1+c4*k2+d4*k3;             if big<abs(y4) return, end
    
42.   k4 = h*f(tj+a4*h,y4);
    
43.   y5 = yj+b5*k1+c5*k2+d5*k3+e5*k4;       if big<abs(y5) return, end
    
44.   k5 = h*f(tj+a5*h,y5);
    
45.   y6 = yj+b6*k1+c6*k2+d6*k3+e6*k4+f6*k5; if big<abs(y6) return, end
    
46.   k6 = h*f(tj+a6*h,y6);
    
47.   err = abs(r1*k1+r3*k3+r4*k4+r5*k5+r6*k6);
    
48.   ynew = yj+n1*k1+n3*k3+n4*k4+n5*k5;
    
49.   if ((err<tol)||(h<2*hmin)),
    
50.     Y(j+1) = ynew;
    
51.     if ((tj+h)>br),
    
52.       T(j+1) = b;
    
53.     else
    
54.       T(j+1) = tj + h;
    
55.     end
    
56.     j = j+1;
    
57.     tj = T(j);
    
58.   end
    
59.   if (err==0),
    
60.     s = 0;
    
61.   else
    
62.     s = 0.84*(tol*h/err)^(0.25);% 最佳步长值
    
63.   end
    
64.   if ((s<0.75)&&(h>2*hmin)), h = h/2; end
    
65.   if ((s>1.50)&&(2*h<hmax)), h = 2*h; end
    
66.   if ((big<abs(Y(j)))||(max1==j)), return, end
    
67. end
    
68. R=[T' Y'];

复制代码

经我测试，结果均准确。后一个程序的精度较高。同步长大致相差一个数量级。

xjzuo

原帖由 bainhome 于 2007-6-9 14:32 发表

                               
登录/注册后可看大图

定步长RK4(自编):
function Y=RungeKutta4(f,xn,y0)
% xn=0:.1:1;

不愧是bainhome , 第二段代码的确处理得很好－－－虽然我以前也用过类似的自适应步长RK法．

[ 本帖最后由 eight 于 2007-6-9 23:47 编辑 ]