二阶常微分方程-龙格库塔法

当刚体受力以后，通常会在介质中运动受到阻力，一般阻力大小正比速度，即kv，或者速度的平方kv^2，方向相反，速度是物体位移-时间微分曲线x(t)的一阶导x’=v，加速度是二阶导x’’=a，假设受力为T，质量为m，这样t时刻，对象受力情况为
f(t)=T-kv=ma
即
T-kx’=mx’’ 化为
x’’=T/m-kx’
这是一个二阶常微分方程
x’’=f(t,x,x’)
初值x(t0)=x0, x’(t0)=x’0

令x’=v，则化为一个一阶常微分方程
v’=f(t,x,v)
四阶龙格库塔法的推导表达式格式如下，其中h为时间步长dt当前步下标为i，下一步为j:
xj=xi + h* vi + (h* h/6.0)(L1+L2+L3);
vj=vi + (h/6.0)(L1+2* L2+2* L3+L4);
其中
hf=h/2.0;
L1=f(ti, xi, vi)
L2=f(ti+hf, xi+hf* vi, vi+hf* L1)
L3=f(ti+hf, xi+hf* vi+hf* hf* L1, vi+hf* L2)
L4=f(ti+h, xi+hf* vi+h* hf* L2, vi+h* L3)