Метод 2-мерного Рунге-Кутты на Mathematica 8 - PullRequest
1 голос
/ 31 декабря 2011

У меня проблема при программировании в 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];

Дело в том, что я пытаюсь реализовать метод Рунге-Кутты.И на самом деле я добился успеха, но только с функцией f[x_], которая имела 1 измерение.Этот код для двух измерений, но он просто не работает, и я не знаю почему.Вот пример кода только с 1-м измерением (обратите внимание, что мне нужно изменить первую строку, чтобы определить функцию, и последнюю строку, когда я вызываю «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];

Чтобы подвести итогКак я реализовал метод Рунге-Кутты для функции с 2 измерениями ??

Если бы вы могли мне помочь, я был бы благодарен.

Заранее спасибо!

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];

Я знаю, в принципе, есть программа математики, которая может решить ЛЮБОЕ уравнение 2-го порядка!Через метод Рунге-Кутта.просто вставьте свою функцию в

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: Задача Рунге-Кутты решает дифференциальные уравнения с проблемами начальных значений, я использовал эту программу в качестве основы для решения проблемы граничных значений, если хотите, просто напишите мне!

1 Ответ

1 голос
/ 09 января 2012

Разве ваш код не реализует то, что уже встроено в Mathematica, а именно, если вы должны использовать опцию

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

для NDSolve?

(Это не значит, чтонет смысла «кататься самостоятельно»: возможно, вы хотите сделать это как учебное упражнение для себя, для студентов или для самого себя.)

...