2
$\begingroup$

I have been asked to reduce the KdV ODE Equation

$d^3/dX^3=(c+6f)df/dX $

to a first order ODE and then solve this ODE by RK4 (Runge-Kutta 4) Numerical Integration.

Analytically, letting c=1 for simplicity, we integrate to get

$f''-3f^2-f=C$

Where $C$ is integration constant. Then Multiply by $f'$ and Integrate to get

$(f')^2/2 -f^3 -f^2/2=Cf +D$

Where $D$ is another integration constant. Now Solve for $f'$ we get

$df/dX=\sqrt(2f^3+f^2+2Cf+2D)$

I have been given initial conditions: $f(0)=-1/2$, $df/dX=0$, and $df^2/dX^2=1/4$

I have been told to:

Use RK4 to integrate numerically the first-order reduction of the ODE from $X = 0$ forward $X_n = n * h, n = 0, 1, . . .N$ and for sufficiently large $X$ (i.e. large $N$). You will need to be plotting $f$ vs $X$ to see if the solution is settling to determine when to stop the integration. Use your judgement as to how small a step size you need to solve this system accurately. If you cannot figure this out from pure thought, experiment with different step sizes and use $δC = |(C(t) − C(t = 0))/C(t = 0)|$ to determine the accuracy. If $δC$ is smaller than $10^{−3}$ for all integration times, then you have acceptable accuracy. Show in a plot your solution for $f(X)$, for $X \in [−20, 20]$.

Here is my Code:

ti = -\[Pi];
tf = \[Pi];
i = 256;
\[CapitalDelta]t = (tf - ti)/i;
X = Range[ti, tf, \[CapitalDelta]t] // N;
F = ConstantArray[0., i + 1];
D1 = ConstantArray[0., i + 1];
D2 = ConstantArray[0., i + 1];
(*Runge Kutta of Order 4*)
F[[1]] = -1/2;
D1[[1]] = 0.;
D2[[1]] = 1/4;

Do[
 k1 = \[CapitalDelta]t*\[Sqrt]Abs[(2 F[[n]]^3 + F[[n]]^2 + 
       2*(D2[[n]] - 3 F[[n]]^2 - F[[n]])* F[[n]] + 
       2*(1/2 D1[[n]]^2 - F[[n]]^3 - 1/2 F[[n]]^2 - D2[[n]] - 
          3 F[[n]]^2 - F[[n]]* F[[n]]))];
 k2 = \[Sqrt]Abs[(2 F[[n]]^3 + F[[n]]^2 + 
       2*(D2[[n]] - 3 F[[n]]^2 - F[[n]])* F[[n]] + 
       2*(1/2 D1[[n]]^2 - F[[n]]^3 - 1/2 F[[n]]^2 - D2[[n]] - 
          3 F[[n]]^2 - F[[n]]* F[[n]]))] + \[CapitalDelta]t*k1/2;
 k3 = \[Sqrt]Abs[(2 F[[n]]^3 + F[[n]]^2 + 
       2*(D2[[n]] - 3 F[[n]]^2 - F[[n]])* F[[n]] + 
       2*(1/2 D1[[n]]^2 - F[[n]]^3 - 1/2 F[[n]]^2 - D2[[n]] - 
          3 F[[n]]^2 - F[[n]]* F[[n]]))] + \[CapitalDelta]t*k2/2;
 k4 = \[Sqrt]Abs[(2 F[[n]]^3 + F[[n]]^2 + 
       2*(D2[[n]] - 3 F[[n]]^2 - F[[n]])* F[[n]] + 
       2*(1/2 D1[[n]]^2 - F[[n]]^3 - 1/2 F[[n]]^2 - D2[[n]] - 
          3 F[[n]]^2 - F[[n]]* F[[n]]))] + \[CapitalDelta]t*k3/2;
 D1[[n + 1]] = D1[[n]] + 1/6*\:f000k1 + 2 k2 + 2 k3 + k4\:f006,
 D2[[n + 1]] = D2[[n]] + 1/6*\:f000k1 + 2 k2 + 2 k3 + k4\:f006,
 F[[n + 1]] = F[[n]] + 1/6*\:f000k1 + 2 k2 + 2 k3 + k4\:f006,
 {n, 1, i}]

So I tried to use 3 vectors with 3 initial values given, and iterate until all the values are filled. but I don't know how to compute the next values of D1 and D2, based on the initial values.

On the other hand, apparently there must be a function with two variable inputs, where it returns a only one F back. I found this online:

Y, T, h = initialize_all(y0, t0, t1, n)
for i in xrange(1, n):
K1 = f(T[i-1], Y[i-1])
tplus = (T[i] + T[i-1]) * .5
K2 = f(tplus, Y[i-1] + .5 * h * K1)
K3 = f(tplus, Y[i-1] + .5 * h * K2)
K4 = f(T[i], Y[i-1] + h * K3)
Y[i] = Y[i-1] + (h / 6.) * (K1 + 2 * K2 + 2 * K3 + K4)
return T, Y

Which is in another language. but clearly my code should follow this pattern. I have tried hard to solve it, but couldn't figure it out. If anyone can help me, I really appreciate their help.

Thanks

$\endgroup$
2
  • 1
    $\begingroup$ Look at your code. There are some funny characters: *\:f000k1" .. Further, change the commas after D1[[n + 1]]..., D2[[n + 1]] =...,F[[n + 1]] =... into semicolons. $\endgroup$ – Daniel Huber Jan 14 at 9:07
  • $\begingroup$ There are so many Runge-Kutta codes on this site it could have its own tag. :) $\endgroup$ – Michael E2 Jan 14 at 18:39
4
$\begingroup$

We should note first that point $X=0$ is a branching point for the first order equation. Therefore we need to make special step before RK4. Using initial condition we calculate, that $C=D=0$, therefore first order equation looks like $$y'=\pm \sqrt {2 y^3+y^2}$$ Code for RK4 we can organize with Module as follows

rk4[f_, h_, x0_, y0_, n_] := 
 Module[{k1, k2, k3, k4}, T[0] = x0; Y[0] = y0;
  Do[k1 = f[T[i - 1], Y[i - 1]]; T[i] = T[i - 1] + h;
   tplus = (T[i] + T[i - 1])*.5;
   k2 = f[tplus, Y[i - 1] + .5*h*k1];
   k3 = f[tplus, Y[i - 1] + .5*h*k2];
   k4 = f[T[i], Y[i - 1] + h*k3];
   Y[i] = Y[i - 1] + (h/6.)*(k1 + 2*k2 + 2*k3 + k4);, {i, 1, n}]; 
  ysol = Table[{T[i], Y[i]}, {i, 0, n}];
  ysol] 

Now we use it in two branches with first step calculation

fp[x_, y_] := Sqrt[2 y^3 + y^2]; h1 = 10^-6; y1 = -1/2 + h1^2/2 *1/4;

yp = rk4[fp, .1, h1, y1, 200];

fm[x_, y_] := -Sqrt[2 y^3 + y^2]; ym = rk4[fm, -.1, -h1, y1, 200];

Finally we plot solution

ListPlot[{ym, yp}, PlotRange -> All, Frame -> True]

Figure 1

To plot function $y(x-t)$ we can use code

ys = Interpolation[Join[ym, yp]];

Plot[Evaluate[Table[ys[x - t], {t, 0, 5, 1}]], {x, -10, 10}, 
 PlotLegends -> Table[t, {t, 0, 5, 1}], Frame -> True]

Figure 2 To compute parameter $C=y''-3 y^2-y$ for numerical solution we should remember that we solve first order ODE, and therefore we need some interpolation of numerical solution to calculate a second derivative. For this problem we use code

uu[h_, h1_, n_] := 
 Module[{hh = h, nn = n}, fp[x_, y_] := Sqrt[2 y^3 + y^2]; 
  y1 = -1/2 + h1^2/2 *1/4;
  yp = rk4[fp, hh, h1, y1, nn]; fm[x_, y_] := -Sqrt[2 y^3 + y^2]; 
  ym = rk4[fm, -hh, -h1, y1, nn]; un = Join[ym, {{0, -1/2}}, yp]; un] 

Let evaluate uu with input data h=h1=.1;n=200 then we have numerical solution as interpolation function

y = Interpolation[uu[.1, .1, 200], InterpolationOrder -> 4]

With this function we can compute C and plot it in a big and small scale as follows (here ch=C but we can't use letter C as a name since it is the default form for the parameter or constant generated in representing the results of various symbolic computations)

ch = y''[x] - 3 y[x]^2 - y[x];

{LogPlot[Abs[ch], {x, -20, 20}, PlotRange -> All, Frame -> True], 
 LogPlot[Abs[ch], {x, -2, 2}, PlotRange -> All, Frame -> True]}

Fugure 3

$\endgroup$
13
  • 1
    $\begingroup$ @AsgharBakerdar See update to my answer. $\endgroup$ – Alex Trounev Jan 14 at 21:03
  • 1
    $\begingroup$ @AsgharBakerdar In general $f(x-t)$ is a traveling wave with waveform $y = f(x)$ traveling to the right (increasing $x$) at unit speed (no matter how $f(x)$ was obtained). $\endgroup$ – Michael E2 Jan 14 at 21:20
  • 1
    $\begingroup$ @MichaelE2 It is wander that NDSolve with Automatic options can't find right numerical solution for this KdV on $(-20,20)$. So we need something very especial like StartingStepSize -> 10^-6, Method -> "ExplicitRungeKutta", WorkingPrecision -> 30, MaxSteps -> Infinity. $\endgroup$ – Alex Trounev Jan 14 at 21:35
  • 1
    $\begingroup$ Well, the constant function $y=-1/2$ is a solution to $y'=\pm \sqrt {2 y^3+y^2}$, so it's hard to blame NDSolve. A standard approach to resolving the singularity at $y=0$ is to differentiate the ODE and use $y'(0)$ and $y''(0)$ to determine the desired branch; but that leads back to where the OP started. -- BTW, my first thought on reading the OP was that the task was to convert KdV to a first-order system in the standard way and apply a vector form of RK4. But you show it's more interesting this way, so I'm not sure what the original task is supposed to be. $\endgroup$ – Michael E2 Jan 14 at 22:02
  • 1
    $\begingroup$ @AsgharBakerdar Since C=y''[x] - 3 y[x]^2 - y[x] and we solve first order equation, we can't use C as a measure of error. We can compute some interpolation for numerical solution and then C. $\endgroup$ – Alex Trounev Jan 14 at 22:31

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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