![](https://csdnimg.cn/release/download_crawler_static/86401249/bg10.jpg)
学院:土木工程与力学学院 姓名:张风顺 学号:320170930401
real(4) k(4,4)!龙格库塔中的 K 矩阵,其 1、2、3、4 行分别指向 x,y,z,w 迭代过程中的 k
值.
f(yy)=u/ak*log(yy/y_0)/(32*d)*c!语句函数
open(1,file="E:\我的作业\计算力学\第三次\四阶龙格库塔\四阶龙格库塔.txt")
y_0=d/30!y0 表示沙粒的起始高度。
do n=10,100,10
write(1,"(a2,i3)")"n=",n
t=0
x=0;y=y_0;z=2.0;w=3.0!初始条件.
call print_answer(t,x,y,z,w)
t=0.6!求解的总时间
t=t/n!时间步长
do i=1,n!求解每个点上的值
!第一层 k 值
k(1,1)=z;k(2,1)=w;k(3,1)=f(y)-c*z/(32*d);k(4,1)=-g-c*w/(32*d)
!第二层 k 值
k(1,2)=z+t/2*k(3,1)
k(2,2)=w+t*k(4,1)/2
k(3,2)=f(y+t/2*k(2,1))-c*(z+t/2*k(3,1))/(32*d)
k(4,2)=-g-c*(w+t/2*k(4,1))/(32*d)
!第三层 k 值
k(1,3)=z+t/2*k(3,2)
k(2,3)=w+t/2*k(4,2)
k(3,3)=f(y+t/2*k(2,2))-c/(32*d)*(z+t/2*k(3,2))
k(4,3)=-g-c/(32*d)*(w+t/2*k(4,2))
!第四层 k 值
k(1,4)=z+t*k(3,3)
k(2,4)=w+t*k(4,3)
k(3,4)=f(y+t*k(2,3))-(z+t*k(3,3))*c/(32*d)
k(4,4)=-g-c/(32*d)*(w+t/2*k(4,3))
!更新 x,y,z,w 的值
x=x+t*(k(1,1)+2*k(1,2)+2*k(1,3)+k(1,4))/6
y=y+t*(k(2,1)+2*k(2,2)+2*k(2,3)+k(2,4))/6
z=z+t*(k(3,1)+2*k(3,2)+2*k(3,3)+k(3,4))/6
w=w+t*(k(4,1)+2*k(4,2)+2*k(4,3)+k(4,4))/6
call print_answer(i*t,x,y,z,w)
end do
end do
contains
subroutine print_answer(t,x,y,z,w)!将 t 时刻的位移和速度打印到指定文件中
real(4) t,x,y,z,w
write(1,*)t,",",x,",",y,",",z,",",w
end subroutine print_answer
end