龙格库塔法求助（万急）

煜宸0922

回复 11 # meiyongyuandeze 的帖子

大侠，再求助一下，我做了边界条件的判定，但是出来的图和原图没有区别，问了下别人，说我的数值精度不对，程序如下，请大侠帮忙指正下错误：
if X(1,j)==0.1
     X(3,j+1)=-0.8*X(3,j);
elseif X(1,j)==-0.1
     X(3,j+1)=-0.8*X(3,j);
end

煜宸0922

回复 6 # meiyongyuandeze 的帖子

大侠，你给的这个程序能不能帮我先改为求下一个值得，不循环。

meiyongyuandeze

回复 17 # 煜宸0922 的帖子

1. clc
   
2. clear
   
3. k=0.001;
   
4. T=0.0;
   
5. A=[0.005 0 0 0]'
   
6. k1=k*feval('fun',T,A);
   
7. k2=k*feval('fun',T+k/2,A+k1/2);
   
8. k3=k*feval('fun',T+k/2,A+k2/2);
   
9. k4=k*feval('fun',T+k,A+k3);
   
10. A_new=A+(k1+2*k2+2*k3+k4)/6

复制代码

比如你知道这一时刻T的变量A的值，求下一时刻T+k时的值A_new，上面的语句就可以实现。

1. A =
   
2.     0.0050
   
3.          0
   
4.          0
   
5.          0
   
6. 
   
7. 
   
8. A_new =
   
9.     0.0050
   
10.     0.0000
    
11.    -0.0009
    
12.     0.0006

复制代码

煜宸0922

回复 18 # meiyongyuandeze 的帖子

十分感谢大侠，我就是在用这个做变步长的龙格库塔的时候老出错，前面有帖子，大侠能不能给指点一二？

煜宸0922

回复 18 # meiyongyuandeze 的帖子

大侠，还有一点时，变步长时有个边界条件abs(X(1,j)-0.1)<=10e-6; X(3,j+1)=-0.8*X(3,j);怎么加在程序里才行呀，我加进去好像老是没有作用，务必请大侠帮助

meiyongyuandeze

回复 20 # 煜宸0922 的帖子

你要加在'fun'这个状态方程中，主程序没调用一次'fun'就会判断你的边界条件！

煜宸0922

回复 21 # meiyongyuandeze 的帖子

可是在fun中，表示碰撞前后速度变化X(3,j+1)=-0.8*X(3,j);没办法表示呀!

meiyongyuandeze

回复 22 # 煜宸0922 的帖子

如果判定出现碰撞，你可以在‘fun’中将其对应的项令为零，然后在主程序中运行X(3,j+1)=-0.8*X(3,j);。

煜宸0922

回复 14 # meiyongyuandeze 的帖子

大侠，麻烦你给我一个完整的变步长龙格库塔的我那个方程组的，我编的老是出问题，怎么想都想不明白！

meiyongyuandeze

我没有现成的程序，不好意思！

煜宸0922

回复 2 # meiyongyuandeze 的帖子

function  y=rk44(x,h,t)
n=1;
k1=h*feval('fun',t(n),x(:,n));
k2=h*feval('fun',t(n)+h/2,x(:,n)+k1/2);
k3=h*feval('fun',t(n)+h/2,x(:,n)+k2/2);
k4=h*feval('fun',t(n)+h,x(:,n)+k3);
y(:,n+1)=x(:,n)+(k1+2*k2+2*k3+k4)/6;
end
大侠，怎么在主程序中调用这个龙格库塔法算的结果进行循环呀？求教!!!

meiyongyuandeze

回复 15 # 煜宸0922 的帖子

给你推荐本书吧，《数值方法（matlab版）》，学校图书馆里应该有

spsmysp

也在做这方面的，学习。

0810064

真是高人楼主走运碰见好人和高人了。

eizhihua602

顶一个。。