Confusion on Runge-Kutta method

Question

Jesse 2013-11-25

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/107488-confusion-on-runge-kutta-method

回答： Meysam Mahooti 2021-5-5

Greetings,

I am trying to implement a 4th order Runge Kutta for the following equations on the interval 1 less than or equal to t less than or equal to 2, with h = 0.01:
   x' = x^-2 + log y + t^2
   y' = e^y-cos(x) + (sin t)x - 1/x^3y^3

with x(2)=-2 and y(2) = 1.

Now, I started to implement it in code is a basic form:

function rungekutta
h = 0.01;
t = 1;
w = 1; 
fprintf('Step 0: t = %12.8f, w = %12.8f\n', t, w);
for i = 1:100
    k1 = h*f(t,w);
    k2 = h*f(t+h/2, w+k1/2);
    k3 = h*f(t+h/2, w+k2/2);
    k4 = h*f(t+h, w+k3);
    w = w + (k1+2*k2+2*k3+k4)/6;
    t = t+h;
    fprintf('Step %d: t = %6.4f, w = %18.15f\n', i, t, w);
end

but the tricky part is that I'm dealing with x,y, and t in both equations, but in searches and books I see, all I see is a function dependent on x (or y) and t with an initial condition, not all three at once.

Any hints or suggestions? I'm a bit stuck in neutral.

Thanks!

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Titus Edelhofer 2013-11-25

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/107488-confusion-on-runge-kutta-method#answer_116267

Hi,

your function f should have two inputs, namely t and X, where X is a two component vector [x, y]. Therefore your w=1 should in fact be w = [-2 1];

Titus

8 个评论
显示 6更早的评论隐藏 6更早的评论

Jesse 2013-12-2

在 MATLAB Online 中打开

Titus, once again thanks for the help. It is now outputting, but not quite correctly.

Once again, my code is now:

function rungekutta
h = 0.01;
t = 1;
w = [-2 1]';
fprintf('Step 0: t = %12.8f, w = %12.8f\n', t, w, w(2));
for i = 1:100
    k1 = h*f(t,w);
    k2 = h*f(t+h/2, w+k1/2);
    k3 = h*f(t+h/2, w+k2/2);
    k4 = h*f(t+h, w+k3);
    w = w + (k1+2*k2+2*k3+k4)/6;
    t = t+h;
    fprintf('Step %d: t = %6.4f, w = %18.15f\n', i, t, w, w(2));
end
function dX = f(t, X)
x = X(1);
y = X(2);
dX(1,1) = x^-2 + log(y) + t^2;
dX(2,1) = exp(y) - cos(x) + sin(t)*x;

and my output is:

>> rungekutta
Step 0: t =   1.00000000, w =  -2.00000000
Step 0: t =   1.00000000, w =   1.00000000
Step 1: t = 1.0100, w = -1.987311093399863
Step 1.014656e+00: t = 1.0147, w = Step 2: t = 1.0200, w = -1.974241485450205
Step 1.029607e+00: t = 1.0296, w = Step 3: t = 1.0300, w = -1.960786707708695
Step 1.044868e+00: t = 1.0449, w = Step 4: t = 1.0400, w = -1.946942103086868
Step 1.060457e+00: t = 1.0605, w = Step 5: t = 1.0500, w = -1.932702816404172
Step 1.076394e+00: t = 1.0764, w = Step 6: t = 1.0600, w = -1.918063783920918
Step 1.092698e+00: t = 1.0927, w = Step 7: t = 1.0700, w = -1.903019721735751
Step 1.109389e+00: t = 1.1094, w = Step 8: t = 1.0800, w = -1.887565112916987
Step 1.126492e+00: t = 1.1265, w = Step 9: t = 1.0900, w = -1.871694193218266
Step 1.144029e+00: t = 1.1440, w = Step 10: t = 1.1000, w = -1.855400935206855
Step 1.162026e+00: t = 1.1620, w = Step 11: t = 1.1100, w = -1.838679030607194
Step 1.180511e+00: t = 1.1805, w = Step 12: t = 1.1200, w = -1.821521870631886
Step 1.199513e+00: t = 1.1995, w = Step 13: t = 1.1300, w = -1.803922524036631
Step 1.219063e+00: t = 1.2191, w = Step 14: t = 1.1400, w = -1.785873712593289
Step 1.239195e+00: t = 1.2392, w = Step 15: t = 1.1500, w = -1.767367783624997
Step 1.259947e+00: t = 1.2599, w = Step 16: t = 1.1600, w = -1.748396679187314
Step 1.281356e+00: t = 1.2814, w = Step 17: t = 1.1700, w = -1.728951901407497
Step 1.303466e+00: t = 1.3035, w = Step 18: t = 1.1800, w = -1.709024473407588
Step 1.326324e+00: t = 1.3263, w = Step 19: t = 1.1900, w = -1.688604895132366
Step 1.349980e+00: t = 1.3500, w = Step 20: t = 1.2000, w = -1.667683093276215
Step 1.374488e+00: t = 1.3745, w = Step 21: t = 1.2100, w = -1.646248364347705
Step 1.399909e+00: t = 1.3999, w = Step 22: t = 1.2200, w = -1.624289309720087
Step 1.426310e+00: t = 1.4263, w = Step 23: t = 1.2300, w = -1.601793761280453
Step 1.453763e+00: t = 1.4538, w = Step 24: t = 1.2400, w = -1.578748695997742
Step 1.482349e+00: t = 1.4823, w = Step 25: t = 1.2500, w = -1.555140137363844
Step 1.512156e+00: t = 1.5122, w = Step 26: t = 1.2600, w = -1.530953041201340
Step 1.543284e+00: t = 1.5433, w = Step 27: t = 1.2700, w = -1.506171162747074
Step 1.575845e+00: t = 1.5758, w = Step 28: t = 1.2800, w = -1.480776901173990
Step 1.609962e+00: t = 1.6100, w = Step 29: t = 1.2900, w = -1.454751116751381
Step 1.645777e+00: t = 1.6458, w = Step 30: t = 1.3000, w = -1.428072914592897
Step 1.683448e+00: t = 1.6834, w = Step 31: t = 1.3100, w = -1.400719387300320
Step 1.723158e+00: t = 1.7232, w = Step 32: t = 1.3200, w = -1.372665306635487
Step 1.765114e+00: t = 1.7651, w = Step 33: t = 1.3300, w = -1.343882751437127
Step 1.809558e+00: t = 1.8096, w = Step 34: t = 1.3400, w = -1.314340655045553
Step 1.856768e+00: t = 1.8568, w = Step 35: t = 1.3500, w = -1.284004250066915
Step 1.907074e+00: t = 1.9071, w = Step 36: t = 1.3600, w = -1.252834380742386
Step 1.960862e+00: t = 1.9609, w = Step 37: t = 1.3700, w = -1.220786642483119
Step 2.018597e+00: t = 2.0186, w = Step 38: t = 1.3800, w = -1.187810292728014
Step 2.080840e+00: t = 2.0808, w = Step 39: t = 1.3900, w = -1.153846854693862
Step 2.148275e+00: t = 2.1483, w = Step 40: t = 1.4000, w = -1.118828301762825
Step 2.221758e+00: t = 2.2218, w = Step 41: t = 1.4100, w = -1.082674658389474
Step 2.302365e+00: t = 2.3024, w = Step 42: t = 1.4200, w = -1.045290771729469
Step 2.391487e+00: t = 2.3915, w = Step 43: t = 1.4300, w = -1.006561875538097
Step 2.490963e+00: t = 2.4910, w = Step 44: t = 1.4400, w = -0.966347344618283
Step 2.603291e+00: t = 2.6033, w = Step 45: t = 1.4500, w = -0.924471646126616
Step 2.731993e+00: t = 2.7320, w = Step 46: t = 1.4600, w = -0.880710769935562
Step 2.882267e+00: t = 2.8823, w = Step 47: t = 1.4700, w = -0.834770995015004
Step 3.062247e+00: t = 3.0622, w = Step 48: t = 1.4800, w = -0.786253802858492
Step 3.285753e+00: t = 3.2858, w = Step 49: t = 1.4900, w = -0.734593453178136
Step 3.579249e+00: t = 3.5792, w = Step 50: t = 1.5000, w = -0.678932918344619
Step 4.004426e+00: t = 4.0044, w = Step 51: t = 1.5100, w = -0.617821940471794
Step 4.776057e+00: t = 4.7761, w = Step 52: t = 1.5200, w = -0.547845348689532
Step 1.276931e+01: t = 12.7693, w = Step 53: t = 1.5300, w =                Inf
Step NaN: t =    NaN, w = Step 54: t = 1.5400, w =                NaN
Step NaN: t =    NaN, w = Step 55: t = 1.5500, w =                NaN
Step NaN: t =    NaN, w = Step 56: t = 1.5600, w =                NaN
Step NaN: t =    NaN, w = Step 57: t = 1.5700, w =                NaN
Step NaN: t =    NaN, w = Step 58: t = 1.5800, w =                NaN
Step NaN: t =    NaN, w = Step 59: t = 1.5900, w =                NaN
Step NaN: t =    NaN, w = Step 60: t = 1.6000, w =                NaN
Step NaN: t =    NaN, w = Step 61: t = 1.6100, w =                NaN
Step NaN: t =    NaN, w = Step 62: t = 1.6200, w =                NaN
Step NaN: t =    NaN, w = Step 63: t = 1.6300, w =                NaN
Step NaN: t =    NaN, w = Step 64: t = 1.6400, w =                NaN
Step NaN: t =    NaN, w = Step 65: t = 1.6500, w =                NaN
Step NaN: t =    NaN, w = Step 66: t = 1.6600, w =                NaN
Step NaN: t =    NaN, w = Step 67: t = 1.6700, w =                NaN
Step NaN: t =    NaN, w = Step 68: t = 1.6800, w =                NaN
Step NaN: t =    NaN, w = Step 69: t = 1.6900, w =                NaN
Step NaN: t =    NaN, w = Step 70: t = 1.7000, w =                NaN
Step NaN: t =    NaN, w = Step 71: t = 1.7100, w =                NaN
Step NaN: t =    NaN, w = Step 72: t = 1.7200, w =                NaN
Step NaN: t =    NaN, w = Step 73: t = 1.7300, w =                NaN
Step NaN: t =    NaN, w = Step 74: t = 1.7400, w =                NaN
Step NaN: t =    NaN, w = Step 75: t = 1.7500, w =                NaN
Step NaN: t =    NaN, w = Step 76: t = 1.7600, w =                NaN
Step NaN: t =    NaN, w = Step 77: t = 1.7700, w =                NaN
Step NaN: t =    NaN, w = Step 78: t = 1.7800, w =                NaN
Step NaN: t =    NaN, w = Step 79: t = 1.7900, w =                NaN
Step NaN: t =    NaN, w = Step 80: t = 1.8000, w =                NaN
Step NaN: t =    NaN, w = Step 81: t = 1.8100, w =                NaN
Step NaN: t =    NaN, w = Step 82: t = 1.8200, w =                NaN
Step NaN: t =    NaN, w = Step 83: t = 1.8300, w =                NaN
Step NaN: t =    NaN, w = Step 84: t = 1.8400, w =                NaN
Step NaN: t =    NaN, w = Step 85: t = 1.8500, w =                NaN
Step NaN: t =    NaN, w = Step 86: t = 1.8600, w =                NaN
Step NaN: t =    NaN, w = Step 87: t = 1.8700, w =                NaN
Step NaN: t =    NaN, w = Step 88: t = 1.8800, w =                NaN
Step NaN: t =    NaN, w = Step 89: t = 1.8900, w =                NaN
Step NaN: t =    NaN, w = Step 90: t = 1.9000, w =                NaN
Step NaN: t =    NaN, w = Step 91: t = 1.9100, w =                NaN
Step NaN: t =    NaN, w = Step 92: t = 1.9200, w =                NaN
Step NaN: t =    NaN, w = Step 93: t = 1.9300, w =                NaN
Step NaN: t =    NaN, w = Step 94: t = 1.9400, w =                NaN
Step NaN: t =    NaN, w = Step 95: t = 1.9500, w =                NaN
Step NaN: t =    NaN, w = Step 96: t = 1.9600, w =                NaN
Step NaN: t =    NaN, w = Step 97: t = 1.9700, w =                NaN
Step NaN: t =    NaN, w = Step 98: t = 1.9800, w =                NaN
Step NaN: t =    NaN, w = Step 99: t = 1.9900, w =                NaN
Step NaN: t =    NaN, w = Step 100: t = 2.0000, w =                NaN
Step NaN: t =    NaN, w = >>