mathematic用龙格库塔法解二阶方程

Question

写一个四阶龙格库塔求解程序
mx''[t]+cx'[t]+kx[t]=-F[t],参数m c k F[t]以及初试条件x[0]=0,x'[0]=0已知，求x[t]
网上我有看到这样的，可以参照。但是我对应解不出来
RungeKutta[func_List, yinit_List, y_List, step_] :=
Module[{k1, k2, k3, k4},
k1 = step N[func /. MapThread[Rule, {y, yinit}]];
k2 = step N[func /. MapThread[Rule, {y, k1/2 + yinit}]];
k3 = step N[func /. MapThread[Rule, {y, k2/2 + yinit}]];
k4 = step N[func /. MapThread[Rule, {y, k3 + yinit}]];
yinit + Total[{k1, 2 k2, 2 k3, k4}]/6]

NestList[RungeKutta[func, #, y, step] &, N[yinit], Round[t/step]]

ListPlot[%, PlotRange -> {{1, 5}, All}]
如有会解方程，可以给我个案例，我可以模仿求解，万分感谢

Answer 1

……所谓龙格库塔法，通俗地说，就是把一个n阶的常微分方程，整理成n个形如 f'(t)=g(t,f(t)) （注意此时右侧不含 f(t) 的导数）的一阶常微分方程组再加以求解的方法。你的方程整理成龙格库塔所需要的形式就是：

x'[t] = y[t]

y'[t] = (-c y[t] - k x[t] - F[t])/m

再考虑到你找到的这段代码本身对格式的要求，就有：

m = 1; c = 1; k = 1; f = 1;
stepsize = 0.01; t = 1;

sol = NestList[
   RungeKutta[{y, (-c y - k x - f)/m}, #, {x, y}, stepsize] &, {0., 0.}, Round[t/stepsize]];

剩下的你自己弄吧……算了，考虑到你的方程里的F[t]涉及了t的具体值，又多了一点点坎（在要使用你找到的这个函数的前提下），再来一个F[t_]=Sin[t]的代码示例吧：

Clear[f, t]
m = 1; c = 1; k = 1;
f[t_] = Sin[t];
stepsize = 0.01; tmin = 0; tmax = 1;

sol = FoldList[
   RungeKutta[{y, (-c y - k x - f[#2])/m}, #, {x, y}, 
     stepsize] &, {0., 0.}, Range[tmin, tmax, stepsize]];