Confusion on Runge-Kutta method

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!

Answer 1

[Titus Edelhofer]

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

[Jesse]

Titus,

Thanks for the response, but there's still more confusion. So, I understand that w should be [-2 1], but a.) does it make a difference that the initial conditions of x and y both occur at 2, and b.) you said the function f should have two inputs.

And this is where I get confused though. See, I have two functions, z' and y' from above, so how would I implement that? I tried adjusting the code to 

function rungekutta(x,y)

   and defining x and y at the command prompt, but I received an error.

   Any suggestions?

   Thanks!

[Titus Edelhofer]

Hi Jesse,

you need to write the function f, e.g. as a subfunction, so add at the bottom:

function dX = f(t, X)
x = X(1);
y = Y(1);
dX(1,1) = ... (that's x')
dX(2,1) = ... (that's y')

Titus

[Jesse]

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!

[Titus Edelhofer]

Hi,

you are nearly there, I guess. Replace in your function f the third line by

y = X(2);

instead of

y = Y(1);

By the way: subfunctions are not much different from functions in other files. The only difference: you can only call them from within the file, so if you try to call it from command window e.g., you will get an error "undefinded function or variable f".

Titus

[Jesse]

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?

[Titus Edelhofer]

Hi,

first, to correct the output, you should change the line to

fprintf('Step %d: t = %6.4f, w = %18.15f, %18.15f\n', i, t, w(1), w(2));

Second, I guess you function f is wrong. At least it's not exactly what you wrote in your original post. I tried using one of MATLAB built-in ODE solvers and they "fail" as well, because the function f "explodes".

Titus

[Jesse]

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

[Steve Chuks]

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]

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