1
$\begingroup$

I am trying to solve y''+4y'+3y=0; y[0]=1, y'[0]=-1 with Runge-Kutta method and I have a problem that its DSolve gives the result:

s = DSolve[{y''[x] + 4 y'[x] + 3 y[x] == 0, y[0] == 1, y'[0] == -1},  y[x], x]

Plot[Evaluate[y[x] /. s], {x, 0, 10}, PlotStyle -> Automatic]

enter image description here

And my code gives a very strange not correct result.

n = 100;
h = 10/n;

y = Table[0, n + 1]; z = Table[0, n + 1];
xxx = Table[0, n + 1];

xxx[[1]] = 0;   y[[1]] = 1;  z[[1]] = -1;

fy[xe_, ye_, ze_] = ze;
fz[x_, yt_, zt_] = -3 yt - 4 zt;


For[i = 1, i <= n, i++,

 Clear[k1z, k2z, k3z, k1y, k2y, k3y];
 yy = Solve[k1y == fy[x + h*1/3, y[[i]] + 1/3*h*k1y, z[[i]]] &&
    k2y == fy[x + h*1, y[[i]] + 1/3*h*k1y + 2/3*h*k2y, z[[i]]] &&
    k3y == fy[x + h*1, y[[i]] + 0 + 0 + 1*h*k3y, z[[i]]], {k1y, k2y, 
    k3y}];


 k1y0 = k1y /. yy[[1]]; k2y0 = k2y /. yy[[1]]; k3y0 = k3y /. yy[[1]];


 zz = Solve[k1z == fy[x + h*1/3, y[[i]], z[[i]] + 1/3*h*k1z] &&
    k2z == fy[x + h*1, y[[i]], z[[i]] + 1/3*h*k1z + 2/3*h*k2z] &&
    k3z == fy[x + h*1, y[[i]], z[[i]] + 0 + 0 + 1*h*k3z], {k1z, k2z, 
    k3z}];

 k1z0 = k1z /. zz[[1]]; k2z0 = k2z /. zz[[1]]; k3z0 = k3z /. zz[[1]];


 y[[i + 1]] = y[[i]] + h (3/4 k1y0 + 3/4 k2y0 - 1/2 k3y0);
 z[[i + 1]] = z[[i]] + h (3/4 k1z0 + 3/4 k2z0 - 1/2 k3z0);

 xxx[[i + 1]] = xxx[[i]] + h; ]

Print[ListPlot[Transpose[{xxx, y}]]]

I have in the result:

enter image description here

So obviously my code is wrong, can anybody help me why?

Maybe somebody can apply new Runge-Kutta or another numerical method to solve it, I would be very happy to see, how it is working?

$\endgroup$
5
  • 1
    $\begingroup$ "[A]nother numerical method to solve it": Try NDSolve[]. $\endgroup$ – Michael E2 May 13 '20 at 18:01
  • $\begingroup$ @Szabolcs unfortunately no, because there is first order of DE, and I need RK for the second order. I think that the problem is somewhere here when I apply RK method for two variables (y, y') $\endgroup$ – Yura Holubeu May 13 '20 at 18:59
  • $\begingroup$ @MichaelE2 no, this is not an option for me $\endgroup$ – Yura Holubeu May 13 '20 at 18:59
  • $\begingroup$ My answer there shows a solution for a 2nd order equation (harmonic oscillator). $\endgroup$ – Szabolcs May 13 '20 at 19:35
  • 1
    $\begingroup$ Then what do you mean by "another numerical method to solve it"? $\endgroup$ – Michael E2 May 13 '20 at 19:41
3
$\begingroup$

This code is original, but it can't be used for nonlinear case. Also it is not looks like RK3. Nevertheless after correcting few typos we have

n = 150;
h = 5/n;

y = Table[0, n + 1]; z = Table[0, n + 1];
xxx = Table[0, n + 1];

xxx[[1]] = 0; y[[1]] = 1; z[[1]] = -1;

fy[xe_, ye_, ze_] := ze;
fz[x_, yt_, zt_] := -3 yt - 4 zt;


Do[
 yy = Solve[
   k1y == fy[xxx[[i]] + h*1/3, y[[i]] + 1/3*h*k1y, z[[i]]] && 
    k2y == fy[xxx[[i]] + h*1, y[[i]] + 1/3*h*k1y + 2/3*h*k2y, 
      z[[i]]] && 
    k3y == fy[xxx[[i]] + h*1, y[[i]] + 0 + 0 + 1*h*k3y, z[[i]]], {k1y,
     k2y, k3y}];
 k1y0 = k1y /. yy[[1]]; k2y0 = k2y /. yy[[1]]; k3y0 = k3y /. yy[[1]];
 zz = Solve[
   k1z == fz[xxx[[i]] + h*1/3, y[[i]], z[[i]] + 1/3*h*k1z] && 
    k2z == fz[xxx[[i]] + h*1, y[[i]], 
      z[[i]] + 1/3*h*k1z + 2/3*h*k2z] && 
    k3z == fz[xxx[[i]] + h*1, y[[i]], z[[i]] + 0 + 0 + 1*h*k3z], {k1z,
     k2z, k3z}];
 k1z0 = k1z /. zz[[1]]; k2z0 = k2z /. zz[[1]]; k3z0 = k3z /. zz[[1]];
 y[[i + 1]] = y[[i]] + h (3/4 k1y0 + 3/4 k2y0 - 1/2 k3y0);
 z[[i + 1]] = z[[i]] + h (3/4 k1z0 + 3/4 k2z0 - 1/2 k3z0);
 xxx[[i + 1]] = xxx[[i]] + h;, {i, 1, n}]

Comparison with exact solution $e^{-x}$

Show[Plot[Exp[-x], {x, 0, 5}, PlotRange -> All], 
 ListPlot[Transpose[{xxx, y}], PlotRange -> All, PlotStyle -> Red]] 

Figure 1

$\endgroup$
12
  • $\begingroup$ woow, it seems you are genius! Can you please tell me where was my mistake in my code? $\endgroup$ – Yura Holubeu May 13 '20 at 20:58
  • $\begingroup$ can you please tell me, if instead of fz[x_, yt_, zt_] := -3 yt - 4 zt; I put in your code fz[x_, yt_, zt_] := 11yt +10 zt; everything falls and solution begans to grow? I just have no idea.... $\endgroup$ – Yura Holubeu May 13 '20 at 21:06
  • $\begingroup$ @YuraHolubeu Actually this is your code, not mine. I just made it working. You use x instead of xxx[[i]] and fy instead of fz. Do you try to implement RK3 schema? $\endgroup$ – Alex Trounev May 13 '20 at 22:20
  • $\begingroup$ yes, it was so stupid that I used x, I couldn't find this for hours... But also I am trying to figure out why I had such big numbers?... because frankly speaking my original problem was with fz[x_, yt_, zt_] := 11yt +10 zt; and RK should give Exp[-x]... I have no idea, why I have such a sharp rise... $\endgroup$ – Yura Holubeu May 13 '20 at 22:46
  • $\begingroup$ and unfortunately if I put fz[x_, yt_, zt_] := x Sqrt[yt]; , then everything begins to count for a really long time, while it should take seconds. I have no idea what is wrong with RK???? Maybe you could know? because probabity that I will figure out it on my own is around zero... $\endgroup$ – Yura Holubeu May 13 '20 at 23:45

Not the answer you're looking for? Browse other questions tagged or ask your own question.