24小时热门版块排行榜    

CyRhmU.jpeg
查看: 2179  |  回复: 17

brqhl_ing

银虫 (小有名气)

[交流] 【求助】四阶龙格-库塔方法!!!!已有3人参与

double precision qx,qy,qz,qx0,qy0,qz0,qx1,qy1,qz1,a,b,c,h,t,d
      double precision qkx1,qkx2,qkx3,qkx4,qky1,qky2,qky3,qky4,qkz1,qkz2
     &,qkz3,qkz4
      double precision px,py,pz,px0,py0,pz0,px1,py1,pz1
      double precision pkx1,pkx2,pkx3,pkx4,pky1,pky2,pky3,pky4,pkz1,pkz2
     &,pkz3,pkz4

        double precision k1,k2 !耦合强度
        double precision  delta, eps  ! 误差精度
      a=0.66d0
      b=0.201d0
      c=0.165d0
        d=1d0/0.407d0
      h=1.0d-3
      px0=0.3d0
      py0=0.5d0
      pz0=1.0d0
        qx0=0.8d0
        qy0=1.5d0
        qz0=0.2d0

        eps=1d-4

        open(10,file='delta.dat')
        k1=0.4
        k2=0.2
!         do k1=0.001d0,1.0d0,0.05d0
!             do k2=0.001d0,1.0d0,0.05d0
                                      
             do 30 t=h,2.0d2,h

               qkx1=h*qx(qy0,qz0,a)
             qky1=h*qy(qx0,qy0,qz0,py0,b,c,k1,k2)
             qkz1=h*qz(qx0,qz0,d)
               pkx1=h*px(py0,pz0,a)
             pky1=h*py(px0,py0,pz0,qy0,b,c,k1,k2)
             pkz1=h*pz(px0,pz0,d)
             qx1=qx0+qkx1/2.0d0
             qy1=qy0+qky1/2.0d0
             qz1=qz0+qkz1/2.0d0
                   px1=px0+pkx1/2.0d0
             py1=py0+pky1/2.0d0
             pz1=pz0+pkz1/2.0d0
        qkx2=h*qx(qy1,qz1,a)
        qky2=h*qy(qx1,qy1,qz1,py1,b,c,k1,k2)
        qkz2=h*qz(qx1,qz1,d)
          pkx2=h*px(py1,pz1,a)
        pky2=h*py(px1,py1,pz1,qy1,b,c,k1,k2)
        pkz2=h*pz(px1,pz1,d)
        qx1=qx0+qkx2/2.0d0
        qy1=qy0+qky2/2.0d0
        qz1=qz0+qkz2/2.0d0
          px1=px0+pkx2/2.0d0
        py1=py0+pky2/2.0d0
        pz1=pz0+pkz2/2.0d0
        qkx3=h*qx(qy1,qz1,a)
        qky3=h*qy(qx1,qy1,qz1,py1,b,c,k1,k2)
        qkz3=h*qz(qx1,qz1,d)
          pkx3=h*px(py1,pz1,a)
        pky3=h*py(px1,py1,pz1,qy1,b,c,k1,k2)
        pkz3=h*pz(px1,pz1,d)
        qx1=qx0+qkx3
        qy1=qy0+qky3
        qz1=qz0+qkz3
        px1=px0+pkx3
        py1=py0+pky3
        pz1=pz0+pkz3
        qkx4=h*qx(qy1,qz1,a)
        qky4=h*qy(qx1,qy1,qz1,py1,b,c,k1,k2)
        qkz4=h*fz(qx1,qz1,d)
          pkx4=h*px(py1,pz1,a)
        pky4=h*py(px1,py1,pz1,qy1,b,c,k1,k2)
        pkz4=h*pz(px1,pz1,d)
        qx1=qx0+(qkx1+2.0*qkx2+2.0*qkx3+qkx4)/6.0d0
        qy1=qy0+(qky1+2.0*qky2+2.0*qky3+qky4)/6.0d0
        qz1=qz0+(qkz1+2.0*qkz2+2.0*qkz3+qkz4)/6.0d0
        px1=px0+(pkx1+2.0*pkx2+2.0*pkx3+pkx4)/6.0d0
        py1=py0+(pky1+2.0*pky2+2.0*pky3+pky4)/6.0d0
        pz1=pz0+(pkz1+2.0*pkz2+2.0*pkz3+pkz4)/6.0d0
           qx0=qx1         
           qy0=qy1            
           qz0=qz1
           px0=px1   
         py0=py1         
           pz0=pz1

大家看看我  我要用四阶龙格-库塔方法解一个六变量的方程组(六个方程) 如何去解  这里我写了一点了  
这个是方程组!求大侠帮忙
回复此楼

» 猜你喜欢

» 本主题相关价值贴推荐,对您同样有帮助:

已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

brqhl_ing

银虫 (小有名气)

nonlinear coupling.obj : error LNK2001: unresolved external symbol _FZ@12
我写的出现了这个错误  是什么原因?
2楼2010-07-24 10:00:38
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

sxf2012

木虫 (正式写手)

★ ★
小木虫(金币+0.5):给个红包,谢谢回帖交流
resonant(金币+1):交流信息费:-) 2010-07-24 11:37:10
qkz4=h*fz(qx1,qz1,d)

fz没有定义吧
3楼2010-07-24 11:21:08
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

maomao1210

金虫 (正式写手)

★ ★
小木虫(金币+0.5):给个红包,谢谢回帖交流
余泽成(金币+1):谢谢参与应助! 2010-07-25 10:01:34
可以参考
Hairer E,etc.
Solving ordinary differential equations
里面有详细讲解
4楼2010-07-24 19:40:02
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

brqhl_ing

银虫 (小有名气)

引用回帖:
Originally posted by sxf2012 at 2010-07-24 11:21:08:
qkz4=h*fz(qx1,qz1,d)

fz没有定义吧

看我后面的函数写的是否有问题  编译能过 但运行不了

double precision function qx(y1,z1,a)
      real*8 y1,z1,a
        qx=y1-a*z1
        return
        end

        double precision function qy(x1,y1,z1,y2,b,c,k1,k2)
       double precision x1,y1,z1,y2,b,c,k1,k2
         qy=-x1+2*b*y1+c*z1-k1*(y1-y2)+k2*(y1**2.0-y2**2.0)
           return
        end

        double precision function qz(x1,z1,d)
        double precision x1,z1,d
        qz=d*(x1-z1**3.0+z1)
        return
        end
      
         double precision function px(y2,z2,a)
        double precision y2,z2,a
        px=y2-a*z2
        return
        end
       double precision function py(x2,y2,z2,y1,b,c,k1,k2)
        double precision x2,y2,z2,y1,b,c,k1,k2
         py=-x2+2*b*y2+c*z2+k1*(y1-y2)-k2*(y1**2.0-y2**2.0)
           return
        end

       double precision function pz(x2,z2,d)
        double precision x2,z2,d
        pz=d*(x2-z2**3.0+z2)
        return
        end
5楼2010-07-25 08:20:11
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

sxf2012

木虫 (正式写手)

★ ★
小木虫(金币+0.5):给个红包,谢谢回帖交流
余泽成(金币+1):谢谢参与应助! 2010-07-26 09:36:17
引用回帖:
Originally posted by brqhl_ing at 2010-07-25 08:20:11:

看我后面的函数写的是否有问题  编译能过 但运行不了

double precision function qx(y1,z1,a)
      real*8 y1,z1,a
        qx=y1-a*z1
        return
        end

        double precision func ...

主程序中,fz 应该改成 pz 吧
6楼2010-07-25 22:29:17
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

brqhl_ing

银虫 (小有名气)

引用回帖:
Originally posted by sxf2012 at 2010-07-25 22:29:17:


主程序中,fz 应该改成 pz 吧

这个问题你一讲过,我就发现了,改了过来。编译能通过。还是不能运行出现两个错误  nonlinear coupling.obj : error LNK2001: unresolved external symbol _Z@8
Debug/nonlinear coupling.exe : fatal error LNK1120: 1 unresolved externals
Error executing link.exe.
7楼2010-07-26 10:49:57
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

brqhl_ing

银虫 (小有名气)

引用回帖:
Originally posted by maomao1210 at 2010-07-24 19:40:02:
可以参考
Hairer E,etc.
Solving ordinary differential equations
里面有详细讲解

这个是一本书吗?现在是程序通不过,那个思想我知道的
8楼2010-07-26 11:37:09
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

sxf2012

木虫 (正式写手)

★ ★
小木虫(金币+0.5):给个红包,谢谢回帖交流
resonant(金币+1):谢谢参与! 2010-08-08 19:55:28
上边贴的程序完整么?我全复制了,编译运行没出现问题
是不是还有其他语句,尤其包括 z(……)的
9楼2010-08-08 19:25:42
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

brqhl_ing

银虫 (小有名气)

引用回帖:
Originally posted by sxf2012 at 2010-08-08 19:25:42:
上边贴的程序完整么?我全复制了,编译运行没出现问题
是不是还有其他语句,尤其包括 z(……)的

double precision qx,qy,qz,qx0,qy0,qz0,qx1,qy1,qz1,a,b,c,h,t,d
      double precision qkx1,qkx2,qkx3,qkx4,qky1,qky2,qky3,qky4,qkz1,qkz2
     &,qkz3,qkz4
      double precision px,py,pz,px0,py0,pz0,px1,py1,pz1
      double precision pkx1,pkx2,pkx3,pkx4,pky1,pky2,pky3,pky4,pkz1,pkz2
     &,pkz3,pkz4

        double precision k1,k2 !耦合强度
        double precision  delta, eps  ! 误差精度
      a=0.66d0
      b=0.201d0
      c=0.165d0
        d=1d0/0.407d0
      h=1.0d-3
      px0=0.3d0
      py0=0.5d0
      pz0=1.0d0
        qx0=0.8d0
        qy0=1.5d0
        qz0=0.2d0

        eps=1d-4

        open(10,file='delta.dat')
!        k1=0.4
!        k2=0.2
         do k1=0.001d0,1.0d0,0.05d0
             do k2=0.001d0,1.0d0,0.05d0
                                      
             do 30 t=h,2.0d2,h

               qkx1=h*qx(qy0,qz0,a)
             qky1=h*qy(qx0,qy0,qz0,py0,b,c,k1,k2)
             qkz1=h*qz(qx0,qz0,d)
               pkx1=h*px(py0,pz0,a)
             pky1=h*py(px0,py0,pz0,qy0,b,c,k1,k2)
             pkz1=h*pz(px0,pz0,d)
             qx1=qx0+qkx1/2.0d0
             qy1=qy0+qky1/2.0d0
             qz1=qz0+qkz1/2.0d0
                   px1=px0+pkx1/2.0d0
             py1=py0+pky1/2.0d0
             pz1=pz0+pkz1/2.0d0
        qkx2=h*qx(qy1,qz1,a)
        qky2=h*qy(qx1,qy1,qz1,py1,b,c,k1,k2)
        qkz2=h*qz(qx1,qz1,d)
          pkx2=h*px(py1,pz1,a)
        pky2=h*py(px1,py1,pz1,qy1,b,c,k1,k2)
        pkz2=h*pz(px1,pz1,d)
        qx1=qx0+qkx2/2.0d0
        qy1=qy0+qky2/2.0d0
        qz1=qz0+qkz2/2.0d0
          px1=px0+pkx2/2.0d0
        py1=py0+pky2/2.0d0
        pz1=pz0+pkz2/2.0d0
        qkx3=h*qx(qy1,qz1,a)
        qky3=h*qy(qx1,qy1,qz1,py1,b,c,k1,k2)
        qkz3=h*qz(qx1,qz1,d)
          pkx3=h*px(py1,pz1,a)
        pky3=h*py(px1,py1,pz1,qy1,b,c,k1,k2)
        pkz3=h*pz(px1,pz1,d)
        qx1=qx0+qkx3
        qy1=qy0+qky3
        qz1=qz0+qkz3
        px1=px0+pkx3
        py1=py0+pky3
        pz1=pz0+pkz3
        qkx4=h*qx(qy1,qz1,a)
        qky4=h*qy(qx1,qy1,qz1,py1,b,c,k1,k2)
        qkz4=h*qz(qx1,qz1,d)
          pkx4=h*px(py1,pz1,a)
        pky4=h*py(px1,py1,pz1,qy1,b,c,k1,k2)
        pkz4=h*pz(px1,pz1,d)
        qx1=qx0+(qkx1+2.0*qkx2+2.0*qkx3+qkx4)/6.0d0
        qy1=qy0+(qky1+2.0*qky2+2.0*qky3+qky4)/6.0d0
        qz1=qz0+(qkz1+2.0*qkz2+2.0*qkz3+qkz4)/6.0d0
        px1=px0+(pkx1+2.0*pkx2+2.0*pkx3+pkx4)/6.0d0
        py1=py0+(pky1+2.0*pky2+2.0*pky3+pky4)/6.0d0
        pz1=pz0+(pkz1+2.0*pkz2+2.0*pkz3+pkz4)/6.0d0
           qx0=qx1         
           qy0=qy1            
           qz0=qz1
           px0=px1   
         py0=py1         
           pz0=pz1

        !响应系统的时间序列写出
         
        delta=((z(1,i)-z(4,i))**2.0+(z(2,i)-z(5,i))**2.0+
     &(z(3,i)-z(6,i))**2.0)**0.5

       
30     continue
      if(delta.le.eps)then
        write(10,'(3f15.6)')k1,k2,0
        write(*,'(3f15.6)')k1,k2,0
        else
        write(10,'(3f15.6)')k1,k2,1
        write(*,'(3f15.6)')k1,k2,1
        endif

        enddo
        enddo
       close(10)
        end

        double precision function qx(y1,z1,a)
      real*8 y1,z1,a
        qx=y1-a*z1
        return
        end

        double precision function qy(x1,y1,z1,y2,b,c,k1,k2)
       double precision x1,y1,z1,y2,b,c,k1,k2
         qy=-x1+2*b*y1+c*z1-k1*(y1-y2)+k2*(y1**2.0-y2**2.0)
           return
        end

        double precision function qz(x1,z1,d)
        double precision x1,z1,d
        qz=d*(x1-z1**3.0+z1)
        return
        end
      
         double precision function px(y2,z2,a)
        double precision y2,z2,a
        px=y2-a*z2
        return
        end
       double precision function py(x2,y2,z2,y1,b,c,k1,k2)
        double precision x2,y2,z2,y1,b,c,k1,k2
         py=-x2+2*b*y2+c*z2+k1*(y1-y2)-k2*(y1**2.0-y2**2.0)
           return
        end

       double precision function pz(x2,z2,d)
        double precision x2,z2,d
        pz=d*(x2-z2**3.0+z2)
        return
        end
10楼2010-08-09 18:33:26
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
相关版块跳转 我要订阅楼主 brqhl_ing 的主题更新
普通表情 高级回复(可上传附件)
信息提示
请填处理意见