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-11-27

在 MATLAB Online 中打开

Titus,

Apologies in the delay, but I think I am getting close. I get this error when I run it:

>> rungekutta
Step 0: t =   1.00000000, w =  -2.00000000
Step 0: t =   1.00000000, w = Undefined function 'Y' for input arguments of type 'double'.
Error in rungekutta>f (line 22)
y = Y(1);
Error in rungekutta (line 11)
    k1 = h*f(t,w);

And the code (with your suggestions I have implemented thanks to your most excellent help) is as follows:

function rungekutta
h = 0.01;
t = 1;  
w   %F1 = h*f(t,x)= [-2 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
function dX = f(t, X)
x = X(1);
y = Y(1);
dX(1,1) = x.^-2 + log(y) + t^2;
dX(2,1) = e.^y - cos(x) +sin(t)*x;

I know MATLAB errors aren't the most helpful, but does the "w" need to become an "x"? I've never worked with subfunctions before unfortunately so this is an awesome opportunity to learn.

Thanks!

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 = >>

It on only appears as though the first there lines of output are correct, and then it turns to nonsense.

Any suggestions?

Jesse 2013-12-10

Titus, I did fix this, and thank you for all of your help!

Steve Chuks 2015-3-6

hi Jesse... could you assist me? I've got a similar work as to solve ODE23s using Runge-Kutta method.

i tried running that code you've got and replaced it with an 'S' function but it exploded. was that what you did?

请先登录，再进行评论。

Answer 2

Meysam Mahooti 2021-5-5

0
链接

此回答的直接链接

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

https://www.mathworks.com/matlabcentral/fileexchange/55430-runge-kutta-4th-order?s_tid=srchtitle

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

请先登录，再进行评论。

Confusion on Runge-Kutta method

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

采纳的回答

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

更多回答（1 个）

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

另请参阅

类别

标签

Community Treasure Hunt

Confusion on Runge-Kutta method

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

采纳的回答

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

更多回答（1 个）

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

另请参阅

类别

标签

Community Treasure Hunt

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

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

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