煜宸0922 发表于 2011-5-11 09:24

回复 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 发表于 2011-5-11 21:42

回复 6 # meiyongyuandeze 的帖子

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

meiyongyuandeze 发表于 2011-5-11 21:56

回复 17 # 煜宸0922 的帖子

clc
clear
k=0.001;
T=0.0;
A='
k1=k*feval('fun',T,A);
k2=k*feval('fun',T+k/2,A+k1/2);
k3=k*feval('fun',T+k/2,A+k2/2);
k4=k*feval('fun',T+k,A+k3);
A_new=A+(k1+2*k2+2*k3+k4)/6比如你知道这一时刻T的变量A的值,求下一时刻T+k时的值A_new,上面的语句就可以实现。A =
    0.0050
         0
         0
         0


A_new =
    0.0050
    0.0000
   -0.0009
    0.0006

煜宸0922 发表于 2011-5-12 08:23

回复 18 # meiyongyuandeze 的帖子

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

煜宸0922 发表于 2011-5-12 10:42

回复 18 # meiyongyuandeze 的帖子

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

meiyongyuandeze 发表于 2011-5-12 11:08

回复 20 # 煜宸0922 的帖子

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

煜宸0922 发表于 2011-5-12 11:35

回复 21 # meiyongyuandeze 的帖子

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

meiyongyuandeze 发表于 2011-5-12 11:54

回复 22 # 煜宸0922 的帖子

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

煜宸0922 发表于 2011-5-12 15:41

回复 14 # meiyongyuandeze 的帖子

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

meiyongyuandeze 发表于 2011-5-12 15:44

我没有现成的程序,不好意思!

煜宸0922 发表于 2011-5-16 10:35

回复 2 # meiyongyuandeze 的帖子

functiony=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 发表于 2011-5-19 10:05

回复 15 # 煜宸0922 的帖子

给你推荐本书吧,《数值方法(matlab版)》,学校图书馆里应该有

spsmysp 发表于 2011-11-17 16:18

也在做这方面的,学习。

0810064 发表于 2011-12-3 18:52

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

eizhihua602 发表于 2012-2-4 08:40

顶一个。。
页: 1 [2]
查看完整版本: 龙格库塔法求助(万急)