急！！！FORTRAN的四阶Runge-kutta方法求解微分方程

annbota

题目：用四阶Runge-kutta方法，求解阻尼振动方程md^2x/dt^2+cdx/dt+kx=0,已知质量m=100,劲度系数k=100,阻尼系数c=20,初速v=0，初始位置x=100,数值计算物体的位置和速度。（用fortran语言)
很急！！！
麻烦各位高手了，谢谢！！！

Seventy721

      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[i],v[i],t[i]);
          em1=h*geval(u[i],v[i],t[i]);
          ek2=h*feval(u[i]+k1/3.0,v[i]+m1/3.0,t[i]+h/3.0);
          em2=h*geval(u[i]+k1/3.0,v[i]+m1/3.0,t[i]+h/3.0);
          ek3=h*feval(u[i]+k1/6.0+k2/6.0,v[i]+m1/6.0+m2/6.0,t[i]+h/3.0);
          em3=h*geval(u[i]+k1/6.0+k2/6.0,v[i]+m1/6.0+m2/6.0,t[i]+h/3.0);
          ek4=h*feval(u[i]+k1/8.0+3.0*k3/8.0,v[i]+m1/8.0+3.0*m3/8.0,t[i]+h/2.0);
          em4=h*geval(u[i]+k1/8.0+3.0*k3/8.0,v[i]+m1/8.0+3.0*m3/8.0,t[i]+h/2.0);
          ek5=h*feval(u[i]+k1/2.0-3.0*k3/2.0+2.0*k4,v[i]+m1/2.0-3.0*m3/2.0+2.0*m4,t[i]+h);
          em5=h*geval(u[i]+k1/2.0-3.0*k3/2.0+2.0*k4,v[i]+m1/2.0-3.0*m3/2.0+2.0*m4,t[i]+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[i+1]=u[i]+(k1+4.0*k4+k5)/6.0;
          v[i+1]=v[i]+(m1+4.0*m4+m5)/6.0;
          t[i+1]=t[i]+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

yangxiaokang

用fortran语言不清楚，不过，用matlab求解很容易的

Seventy721

给我你的邮箱我给你发代码。可以发私人消息给我。

lavi

楼上是北师大二附中的么？

凌绝顶

Seventy721 发表于 2010-4-2 13:41

                               
登录/注册后可看大图

给我你的邮箱我给你发代码。可以发私人消息给我。

有好东西贡献出来嘛，还私聊什么啊

Rainyboy

回复 5 # 凌绝顶 的帖子

我没有针对谁，只是想随便说说。
怎么说呢，自己的东西再不完善也是自己的东西，不大喜欢被别人随便用，也不喜欢分享；而别人的好东西呢，可以通过“求算法”，“求调试”等等来简单地得到，这都还是好的了，很多大公司的软件是“求”不到的，怎么办呢，用盗版呗，而且盗版得没有一点愧疚。
很多人，毫不隐晦的说，包括我在内，大概是这么一种心理吧。

赤血冰霜

yangxiaokang 发表于 2010-3-8 21:17

                               
登录/注册后可看大图

用fortran语言不清楚，不过，用matlab求解很容易的

那请问，知道一个含有未知数的行列式值为0，怎么用MATLAB求这个未知数啊？

Rainyboy

回复 7 # 赤血冰霜 的帖子

实际上就是高次方程的求根，你遇到的是什么问题？求特征值么？

赤血冰霜

回复 8 # Rainyboy 的帖子

对，就是求特征值的，

Rainyboy

回复 9 # 赤血冰霜 的帖子

在MATLAB版搜一搜就能找到一大堆啊……用eig函数就行

赤血冰霜

回复 10 # Rainyboy 的帖子

对MATLAB态不了解了，您能不能推荐本书啊，另外，如果知道一个行列式的值，且该行列式含有未知数，怎么求该未知数啊？谢谢了

Rainyboy

回复 11 # 赤血冰霜 的帖子

我也不怎么用MATLAB，不敢乱推荐，呵呵。

含未知数的行列式，其实本身就是一个高次方程，所有可以用在高次方程求根的算法都可以用在这里，不明白你纠结在什么地方？把你的“含未知数的行列式”给出来看看？

赤血冰霜

回复 12 # Rainyboy 的帖子

比如，DET(A)=0,A=[]（5行5列）,其中含一未知数，怎么在MATLAB中实现未知数的求解？

attempter

回复 3 # Seventy721 的帖子

这个我也很需要 能发给我吗 876698559@qq.com