Mathematica 8 上的二维龙格-库塔方法

Question

我在 Mathematica 8 中编程时遇到问题，这是我的代码：

f[t_, y_] := {y, y};

RungeKutta3[a_, b_, Alpha_, n_, f_] := 
  Module[{h, j, k1, k2, k3}, 
    h = (b - a)/n; 
    Y = T = Table[0, {100 + 1}]; 
    Y[[1]] = Alpha; 
    T[[1]] = a; 
    For[j = 1, j <= n, ++j, 
      k1 = f[T[[j]], Y[[j]]]; 
      k2 = f[T[[j]] + h/2, Y[[j]] + k1*h/2]; 
      k3 = f[T[[j]] + h, Y[[j]] + (-k1 + 2 k2)h]; 
      Y[[j + 1]] = Y[[j]] + h/6(k1 + 4 k2 + k3); 
      (* Print[j, "----->", Y[[j]]];*) 
      T[[j + 1]] = T[[j]] + h;
   ];]; 

RungeKutta3[0., 1., {300., 500}, 2, f];

问题是，我正在尝试实现 Runge-Kutta 方法。我实际上是成功的，但只是使用了具有 1 维的函数 f[x_]。该代码适用于二维，但它根本不起作用，我不知道为什么。这是仅具有一维的代码的示例（请注意，当我调用“RungeKutta3”时，我必须更改第一行来定义函数和最后一行）。

f[t_, y_] := y;

RungeKutta3[a_, b_, Alpha_, n_, f_] := 
  Module[{h, j, k1, k2, k3}, 
    h = (b - a)/n;  
    Y = T = Table[0, {100 + 1}];  
    Y[[1]] = Alpha; 
    T[[1]] = a;  
    For[j = 1, j <= n, ++j,   
      k1 = f[T[[j]], Y[[j]]];   
      k2 = f[T[[j]] + h/2, Y[[j]] + k1*h/2];   
      k3 = f[T[[j]] + h, Y[[j]] + (-k1 + 2 k2)*h];   
      Y[[j + 1]] = Y[[j]] + h/6*(k1 + 4 k2 + k3);   
      (* Print[j, "----->", Y[[j]]];*)     
      T[[j + 1]] = T[[j]] + h;
  ];]; 

RungeKutta3[0., 1., 300., 100, f];

综上所述，如何实现二维函数的龙格-库塔方法？

如果你能帮助我，我将不胜感激。

提前致谢！

PS：龙格-库塔方法是3阶

----------------------

问题解决了！检查代码，如果有人需要帮助，尽管问！

f[t_, y1_, y2_] := 3 t*y2 + Log[y1] + 4 y1 - 2 t^2 * y1 - Log[t^2 + 1] - t^2;
F[t_, {y1_, y2_}] := {y2, f[t, y1, y2]}; 
RungeKutta3[a_, b_, [Alpha]_, n_, f_] :=
 Module[{h, j, k1, k2, k3, Y, T, R},
  h = (b - a)/n;
  Y = T = Table[0, {n + 1}];
  Y[[1]] = [Alpha]; T[[1]] = a;
  For[j = 1, j <= n, ++j,
   k1 = f[T[[j]], Y[[j]]];
   k2 = f[T[[j]] + h/2, Y[[j]] + k1*h/2];
   k3 = f[T[[j]] + h, Y[[j]] + (-k1 + 2 k2)*h];
   Y[[j + 1]] = Y[[j]] + h/6*(k1 + 4 k2 + k3);
   T[[j + 1]] = T[[j]] + h;
   ];
  R = Table[0, {n + 1}]; 
  For[j = 1, j <= n + 1, j++, R[[j]] = Y[[j]][[1]]];
  Print[ListPlot[Transpose[{T, R}]]]
  ];

RungeKutta3[0., 1, {1., 0.}, 1000, F];

我知道基本上有一个数学程序可以解决任何二阶方程！通过龙格-库塔法。只需将函数插入其中

f[t_, y1_, y2_]:= [Insert your function here]

t 是独立值，y1 是函数本身 y(t)，y2 是 y'(t)。 通过以下方式调用函数：

RungeKutta3[a, b, [Alpha], n, F];

其中 a 是初始“t”值，b 最终“t”值，[Alpha] 初始值函数和一阶导数（以 {y1(a),y2(a0)} 形式给出），n 是您想要表示的等距点数。 F 是您必须插入的函数，尽管您为 f 提供了函数。

如有任何问题，请随时提出！ PS：龙格-库塔问题解决了带有初始值问题的微分方程，我用这个程序作为解决边界值问题的基础，如果你想要它，请给我发短信！

Answer 1

您的代码是否只实现了 Mathematica 中已内置的内容，即，如果您要使用

Method -> {"ExplicitRungeKutta", "DifferenceOrder" -> 3}

NDSolve 选项？

（这并不是说“自己动手”没有价值：也许您想将其作为自己或学生或作为学生自己的学习练习。）