急！！！FORTRAN的四阶Runge-kutta方法求解微分方程
题目：用四阶Runge-kutta方法，求解阻尼振动方程md^2x/dt^2+cdx/dt+kx=0,已知质量m=100,劲度系数k=100,阻尼系数c=20,初速v=0，初始位置x=100,数值计算物体的位置和速度。（用fortran语言)很急！！！
麻烦各位高手了，谢谢！！！ x0=0
y0=1
xf=1
nn=100
h=(xf-x0)/nn
x=x0
y=y0
c
100 continue
call rk5(x,y,xx,yy)
c
if(x.lt. xf)
x=x+h
go to 100
end if
stop
end
c
c...This subroutine solves nonlienar ODE using the fifth order Rounge-Kutta Method
c suggested bby Butcher - J. Australian Math. Soc., Vol. 4, 1964, o.179.
c
subroutine rk5 (xi,yi,h,yy)
ek1=h*feval(u,v,t);
em1=h*geval(u,v,t);
ek2=h*feval(u+k1/3.0,v+m1/3.0,t+h/3.0);
em2=h*geval(u+k1/3.0,v+m1/3.0,t+h/3.0);
ek3=h*feval(u+k1/6.0+k2/6.0,v+m1/6.0+m2/6.0,t+h/3.0);
em3=h*geval(u+k1/6.0+k2/6.0,v+m1/6.0+m2/6.0,t+h/3.0);
ek4=h*feval(u+k1/8.0+3.0*k3/8.0,v+m1/8.0+3.0*m3/8.0,t+h/2.0);
em4=h*geval(u+k1/8.0+3.0*k3/8.0,v+m1/8.0+3.0*m3/8.0,t+h/2.0);
ek5=h*feval(u+k1/2.0-3.0*k3/2.0+2.0*k4,v+m1/2.0-3.0*m3/2.0+2.0*m4,t+h);
em5=h*geval(u+k1/2.0-3.0*k3/2.0+2.0*k4,v+m1/2.0-3.0*m3/2.0+2.0*m4,t+h);
E1=(2.0*k1-9.0*k3+8.0*k4-k5)/30.0;
E2=(2.0*m1-9.0*m3+8.0*m4-m5)/30.0;
u=u+(k1+4.0*k4+k5)/6.0;
v=v+(m1+4.0*m4+m5)/6.0;
t=t+h;
c
x1=xi
y1=yi
call func(x1,y1,f1)
ek1=h*f1
x2=xi+h/4
y2=yi+ek1/4
call func(x2,y2,f2)
ek2=h*f2
c
x3=xi+h/4
y3=yi+ek1/8+ek2/8
call func(x3,y3,f3)
ek3=h*f2
c
x4=xi+h/2
y4=yi-ek2/2+ek3
call func(x4,y4,f4)
ek4=h*f4
c
x5=xi+3.*h/4.
y5=yi+3.*ek1/16.+9.*ek4/16.
call func(x5,y5,f5)
ek5=h*f5
c
x6=xi+h
y6=yi-3.*ek1/7.+2.*ek2/7.+12.*ek3/7.-12.*ek4/7.+8.*ek6/7.
call func(x6,y6,f6)
ek6=h*f6
c
yy=yi+(7.*ek1+32.*ek3+12.*ek4+32.*ek5+7.*ek6)
return
end
subroutine func(x,y,f,g)
c
c...This subroutine claculates the value of f(x,y)
c
f=x
return
end
回复 楼主 annbota 的帖子
用fortran语言不清楚，不过，用matlab求解很容易的 给我你的邮箱我给你发代码。可以发私人消息给我。 楼上是北师大二附中的么？ Seventy721 发表于 2010-4-2 13:41 static/image/common/back.gif给我你的邮箱我给你发代码。可以发私人消息给我。
有好东西贡献出来嘛，还私聊什么啊 回复 5 # 凌绝顶 的帖子
我没有针对谁，只是想随便说说。
怎么说呢，自己的东西再不完善也是自己的东西，不大喜欢被别人随便用，也不喜欢分享；而别人的好东西呢，可以通过“求算法”，“求调试”等等来简单地得到，这都还是好的了，很多大公司的软件是“求”不到的，怎么办呢，用盗版呗，而且盗版得没有一点愧疚。
很多人，毫不隐晦的说，包括我在内，大概是这么一种心理吧。 yangxiaokang 发表于 2010-3-8 21:17 static/image/common/back.gif
用fortran语言不清楚，不过，用matlab求解很容易的
那请问，知道一个含有未知数的行列式值为0，怎么用MATLAB求这个未知数啊？ 回复 7 # 赤血冰霜 的帖子
实际上就是高次方程的求根，你遇到的是什么问题？求特征值么？ 回复 8 # Rainyboy 的帖子
对，就是求特征值的， 回复 9 # 赤血冰霜 的帖子
在MATLAB版搜一搜就能找到一大堆啊……用eig函数就行 回复 10 # Rainyboy 的帖子
对MATLAB态不了解了，您能不能推荐本书啊，另外，如果知道一个行列式的值，且该行列式含有未知数，怎么求该未知数啊？谢谢了 回复 11 # 赤血冰霜 的帖子
我也不怎么用MATLAB，不敢乱推荐，呵呵。
含未知数的行列式，其实本身就是一个高次方程，所有可以用在高次方程求根的算法都可以用在这里，不明白你纠结在什么地方？把你的“含未知数的行列式”给出来看看？ 回复 12 # Rainyboy 的帖子
比如，DET(A)=0,A=[]（5行5列）,其中含一未知数，怎么在MATLAB中实现未知数的求解？ 回复 3 # Seventy721 的帖子
这个我也很需要 能发给我吗 876698559@qq.com
页:
[1]
2