| 查看: 2323 | 回复: 17 | |||
[交流]
【求助】四阶龙格-库塔方法！！！！ 已有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 大家看看我 我要用四阶龙格-库塔方法解一个六变量的方程组（六个方程） 如何去解 这里我写了一点了  这个是方程组！求大侠帮忙 |
回复此楼» 猜你喜欢
伙伴们，祝我生日快乐吧
已经有25人回复
-
311求调剂
已经有5人回复
-
0703一志愿211 285分求调剂
已经有5人回复
-
070300化学学硕求调剂
已经有6人回复
-
283求调剂
已经有14人回复
-
中科院材料273求调剂
已经有4人回复
-
085600调剂
已经有6人回复
-
0703化学调剂，求各位老师收留
已经有8人回复
-
070303一志愿西北大学学硕310找调剂
已经有8人回复
-
085600材料与化工 求调剂
已经有14人回复
» 本主题相关价值贴推荐，对您同样有帮助:
- 隐式欧拉法求解一阶常微分方程
已经有7人回复
- LABVIEW四阶龙格库塔法求解一阶微分方程组
已经有12人回复
- 【求助】注册化工工程师的报考条件，报考时间以及考试科目都是什么？
已经有26人回复
- 【求助】Matlab中利用四阶龙格-库塔法求解微分方程！！！！
已经有9人回复
- 【讨论】【【四阶龙格-库塔法求解微分方程】】
已经有8人回复
- 注册化工工程师
已经有10人回复
- 【分享】暑假闲来无事，与虫虫共同进步
已经有23人回复
- 【求助】哪位高手帮我翻译下这段英文算法的说明【已解决】
已经有5人回复
- 【求助】使用Matlab拟合反应动力学方程问题
已经有7人回复
2楼2010-07-24 10:00:38
3楼2010-07-24 11:21:08
maomao1210
金虫 (正式写手)
- 程序强帖: 5
- 应助: 2 (幼儿园)
- 金币: 1432.8
- 散金: 242
- 红花: 16
- 沙发: 1
- 帖子: 991
- 在线: 442.1小时
- 虫号: 253215
- 注册: 2006-05-20
- 性别: MM
- 专业: 考古理论
4楼2010-07-24 19:40:02
|
看我后面的函数写的是否有问题 编译能过 但运行不了 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
6楼2010-07-25 22:29:17
7楼2010-07-26 10:49:57
8楼2010-07-26 11:37:09
9楼2010-08-08 19:25:42
|
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
