Runge Kutta method for a system

Questions
,
1

I want to solve x’=y,y’=6x-y+G,x(0)=1,y(0)=1 with Runge-Kutta Method. Here G lies in [0,1]. I want to get the plot of G versus x. is that possible? May any one help me with that? Here is my code.

ClearAll["Global`*"];    
a = 0; b = 1;
x[0] = 1; y[0] = 1;
h = 0.1; n = (b - a)/h;
Do[t[i] = 0.0 + (i)*h, {i, 0, n}]
Do[G[j] = 0.0 + (j)*h, {j, 0, n}]
f[t_, G_,x_, y_] = y; g[t_,G_, x_, y_] = 6*x - y+G;
Do[S1 = f[t[i], G[i], x[i], y[i]]; P1 = g[t[i], G[i], x[i], y[i]]; 
S2 = f[t[i] + h, G[i] + h, x[i] + h*S1, y[i] + h*P1]; 
P2 = g[t[i] + h, G[i] + h, x[i] + h*S1, y[i] + h*P1]; 
S3 = f[t[i] + h/2, G[i] + h/2, x[i] + h*S2/2, y[i] + h*P2/2]; 
P3 = g[t[i] + h/2, G[i] + h/2, x[i] + h*S2/2, y[i] + h*P2/2];
S4 = f[t[i] + h, G[i] + h, x[i] + h*S3, y[i] + h*P3];
P4 = g[t[i] + h, G[i] + h, x[i] + h*S3, y[i] + h*P3];
x[i + 1] = x[i] + (S1 + 2 S2 + 2 S3 + S4)*h/6; 
y[i + 1] = y[i] + (P1 + 2 P2 + 2 P3 + P4)*h/6, {i, 0, k}]
Table[{G[j], x[j]}, {j, 0, k}];
2

What do you mean by “plot of G versus x”? Isn’t G constant?

And what about this?

eq1 = x'[t] == y[t];
eq2 = y'[t] == 6 x[t] - y[t] + g;
DSolve[{eq1, eq2, x[0] == 0, y[0] == 1}, {x, y}, t]
3

Thank you for your reply.G is constant.but let me explain in detail:
Let us assume G=0. one gets the value of x(t) versus this G right?
Again, let us assume G=0.1.one gets the value of x(t) versus this G right?
Again, let us assume G=0.2.one gets the value of x(t) versus this G right?
and so on.
Now, is it logical to examine the behaviolur x(t) versus these G’s?
I mean: is it possible get the plot of G versus x(t)?
Do you think it is well posed problem?Thanks.

4

Hello Vedat, I am afraid I do not understand what you mean. What do you mean by “one gets the value of x[t] versus G” ? x[t] is a function and not a value! x[ .85 ] for example is a value, and this value certainely depends on G. But for a given t x[ t ] is a linear function of G (look for example at fx[.1, g] ), so you get a straight line. The same is true for y,
I attach some code which shows how the functions x[t] and y[t] ( and not a value ) change with G (in my code g)

eq1 = x'[t] == y[t];
eq2 = y'[t] == 6 x[t] - y[t] + g;
sol = DSolve[{eq1, eq2, x[0] == 0, y[0] == 1}, {x, y}, t];
fx[t_, g_] := Evaluate[(x /. sol[[1, 1]])[[2]]]
fy[t_, g_] := Evaluate[(y /. sol[[1, 2]])[[2]]]
(* Show x[t] depending on g  *)
Manipulate[
 Plot[fx[tt, g], {tt, 0, 2}, PlotRange -> {0, 6}],
 {g, 0, 1}]
(* Show both x[t] (red) y[t] (blue)  depending on g  *)
Manipulate[
 Plot[{fx[t, g], fy[t, g]}, {t, 0, 2},
  PlotRange -> {0, 6},
  PlotStyle -> {Red, Blue}],
 {g, 0, 1}]
(* Show the vector { x[t], y[t] } depending on g  *)
Manipulate[
 ParametricPlot[{fx[t, g], fy[t, g]}, {t, 0, 2},
  PlotRange -> {{0, 2}, {0, 4}}],
 {g, 0, 1}]
5

Hello Hans, thank you for your invaluable help,again. I had to set up this simple problem. Because actually my goal was to solve the problem in the attached article. I hope you have time for this. The authors solve the system (4.1) on page 11 with the fourth order Runge-Kutta method.
Here, except x and y, all other parameters are constant.After calculating x(t) and y(t), we can use (4.5) with the help of (4.2). But for this, we have to find out x(t) and y(t). Please see figure 4 now. g versus Q(s). is something like this possible? that’s why I had to impose the simple problem I sent you
Am I not right?

DM.pdf (984 KB)

6

Hello Vedat,

I had a look at the paper. Equation 4.1 for me reads (please check it)

eq11a = x'[t] == y[t]
eq11b = lam (x - lam x^3 + lam^2 x^4) (x'[t])^2 - 
alf gam (1 + lan x^2) x'[t] - w0^2 x - del x^3 - xi x^5 + 
gam (1 + lam x^2) ( f Cos[w t] + g cos[ww t]) /. {x'[t] -> y[t],  x -> x[t]}

Then there are values given like (I omit the beta, it doesn’t appear in the equation)

val = {alf -> .2, w0 -> I, ww -> 9.842, w -> .5, f -> .05}

but inserting these leaves several items open, so a numerical integration ist not possible

eq2 /. val

What about xi, gam and so on?

7

Hello again Hans;
Please see the attached files for my parameters.we can use them to get for example figure 4 c.&By the way, in your eq11b lam^2 x^4 should be replaced with lam^2 x^5Wolfram Notebook

parameters.pdf (117 KB)

8

Hi Vedat,

here is a code which gives you some values of Qs and Qc, which of course depend on g.
Please have a look about typos and if the results are reasonable.

I still don’t know what you finally want to do, x[t] is a function and not a value, but I admit I didn’t read the paper, only the part around equation ( 4.1 )

eq11a = x'[t] == y[t];
eq11b = y'[t] == 
    lam (x - lam x^3 + lam^2 x^5) (x'[t])^2 - 
     alf gam (1 + lam x^2) x'[t] - w0^2 x - del x^3 - xi x^5 + 
     gam (1 + lam x^2) (f Cos[w t] + g Cos[ww t]) /. {x'[t] -> y[t], 
    x -> x[t]};

val = {alf -> .2, w0 -> I, ww -> 9.842, w -> .5, f -> .05, lam -> 0, 
   beta -> 1, m0 -> 1.5, gam -> 1/m0, del -> beta gam - lam w0^2, 
   xi -> beta gam lam +( lam  w0)^2};

eq11b2 = eq11b //. val;

sol[gg_?NumericQ, t1_?NumericQ] := 
 Module[{}, eq11b3 = eq11b2 /. g -> gg;
  NDSolve[{eq11a, eq11b3, x[0] == 1, y[0] == 1}, {x, y}, {t, 0, t1}]]

nT = 10;
lsg = sol[.3, nT];

fx = x /. lsg[[1, 1]]
fy = y /. lsg[[1, 2]]

qQs = NIntegrate[fx[t] (Sin[w t] /. val), {t, 0, nT}] 2/nT
qQc = NIntegrate[fx[t] (Cos[w t] /. val), {t, 0, nT}] 2/nT
9

Now first I should calculate Qnum given in Eq.(4.5) and then get the plot of g versus Qnum
for this m0 given in Fig.4 c color in green interval in [0,200] for g.

10

Hmmmm ?

Ok, add to the foregoing code this (but there seems to difficulties during the integration, I didn’t have a closer look at this) It doesn’t look similar to fig. 4c.

qQnum[g_] := Module[{},
  nT = 10;
  lsg = sol[g, nT];
  fx = x /. lsg[[1, 1]];
  fy = y /. lsg[[1, 2]];
  qQs = NIntegrate[fx[t] (Sin[w t] /. val), {t, 0, nT}] 2/nT;
  qQs = NIntegrate[fx[t] (Sin[w t] /. val), {t, 0, nT}] 2/nT;
  qQc = NIntegrate[fx[t] (Cos[w t] /. val), {t, 0, nT}] 2/nT;
  qnum = Sqrt[qQs^2 + qQc^2]]

and

tt = Table[{qQnum[g], g}, {g, 0, 200}];
ListLinePlot[tt]
ListLinePlot[Reverse /@ tt]
11

I think

Table[{qQnum[g], g}, {g, 0, 200}]

should be replaced with

Table[{g,qQnum[g]}, {g, 0, 200}]

Right?

12

Hmmmm ?

Ok, add to the foregoing code this (but there seems to difficulties during the integration, I didn’t have a closer look at this) It doesn’t look similar to fig. 4c.

qQnum[g_] := Module[{},
  nT = 10;
  lsg = sol[g, nT];
  fx = x /. lsg[[1, 1]];
  fy = y /. lsg[[1, 2]];
  qQs = NIntegrate[fx[t] (Sin[w t] /. val), {t, 0, nT}] 2/nT;
  qQs = NIntegrate[fx[t] (Sin[w t] /. val), {t, 0, nT}] 2/nT;
  qQc = NIntegrate[fx[t] (Cos[w t] /. val), {t, 0, nT}] 2/nT;
  qnum = Sqrt[qQs^2 + qQc^2]]

and

tt = Table[{qQnum[g], g}, {g, 0, 200}];
ListLinePlot[tt]
ListLinePlot[Reverse /@ tt]
13

Reverse@tt interchanges q and qQnum. So you get both plots,

14

Hello Hans, According to the paper, n=1,2,3,…,T=2Pi/w and Ts=nT.
I have added them to the code. But I can not run it.I need your help.

Open in Wolfram Cloud ↗

Bel1.pdf (207 KB)

15

The code works fine for me

tt = Table[{g, qQnum[g]}, {g, 0, Ts}]
(*
{{0, 0.303553}, {1, 0.299951}, {2, 0.296613}, {3, 0.293514},
 {4, 0.290635}, {5, 0.287959}, {6, 0.285475}, {7, 0.283173},
 {8, 0.281044}, {9, 0.279083}, {10, 0.277284}, {11, 0.275644},
 {12, 0.27416}}
*)

ListLinePlot[tt]

10779plot.png
10779plot.png360×220 6.57 KB

16

Helo hans, did you add these n = 1; w = 0.5; T = (2N[Pi])/w; Ts = nT;

17

I can not get the same result Could you please as a single code again?

18

I just copied and pasted all of the code you provided above in a new notebook and ran it against a fresh kernel. Generating the tt table takes ~45 sec. to finish on my laptop running “12.2.0 for Mac OS X x86 (64-bit) (December 12, 2020)”

19

Dear Hans,

Now it is ok. I fixed it. As a last thing, was my question logical?
I mean the plot of G versus qQnum is logical? Right?

20

Hello Vedat,

I am afraid I cannot comment on that. I didn’t read the whole paper you attached some posts above. What do you mean by “Is it logical to plot quantity1 against quantity2” ? Of course it is. The question is, what do you want to show.

I still don’t know what you essentially want to do. Reproduce some figures of the paper? Doing some new things? Please excuse me, but I have no idea what you wanted to do. Rohit says, the code (which version of code exactly) works on his system, and he gets a reasonable plot of qQnum against g. If it this what you wanted I dare say, everything is ok. Or isn’t it?

@Rohit: Did you not encounter error messages during the numerical integration? I got some errors concerning convergence of the integrals.