Numerical solution of KdV equation

Question

I am dealing with the following Cauchy problem for the KdV equation:

w'''[t] + 6 w[t] w'[t] - w'[t] = HeavisideTheta[t]
w[0] = w'[0] = w''[0] = 0

My problem is to minimize the distance between w[t] and wG[t], where

wG[t] = s0 Integrate[G[t1] HeavisideTheta[t - t1],{t, 0, t}]

Here G is the as the general solution of the following Cauchy problem:

G'''[t] + 6 G[t] G'[t] - G'[t] = 0,
G[0] = G'[0] = 0, G''[0] = s.

Here s is a real number.

Thus, I aim to solve

NMinimize[MaxValue[Abs[w[t] - wG[t]]],{s, s0}]

with respect to $s$ and $s_0$ when 0 <= t <= 5.

Note that I do not want to reduce the order of the KdV equation as it is usually done by integrating it with respect to $t$ and neglecting the integration constant.

My attempt. First, I get the numerical solution of the original KdV:

theta[t_] := 0.5 (1 + Tanh[100 t])
f[t_] := theta[t]
solw = NDSolve[{w'''[t] + 6 w[t] w'[t] - w'[t] == f[t], w[0] == 0, 
    w'[0] == 0, w''[0] == 0}, w, {t, 0, 20}];
wsol[t_] := Evaluate[w[t] /. solw]

Then, I get the numerical solution the equation for the Green's function:

solG = ParametricNDSolve[{G'''[t] + 6 G[t] G'[t] - G'[t] == 0, 
    G[0] == 0, G'[0] == 0, G''[0] == s}, G, {t, 0, 20}, {s}, 
   MaxSteps -> Infinity, Method -> {"ParametricCaching" -> None}];
GGreen[t_, s_] := Evaluate[G[s][t] /. solG]

And finally, I compute $\tilde{w}$ by Simpsons's rule:

En = 50;
T[n_] := t/En n
GsolSimp[t_, s_, s0_] := 
 s0 t/(3 En) (GGreen[T[0], s] f[t - T[0]] + 2 Sum[GGreen[T[2  n], s] f[
        t - T[2  n]],{n, 1, En/2 - 1}] + 4 Sum[GGreen[T[2  n - 1], s] f[
        t - T[2  n - 1]],{n, 1, En/2}] + GGreen[T[En], s] f[t - T[En]])

When I run the command

NMinimize[{Max[
   Table[Abs[wsol[t] - GsolSimp[t, s, s0]], {t, 0, 5, 0.1}]], 
  0 <= s <= 1, -10 <= s0 <= 10}, {s, s0}, 
 Method -> {"RandomSearch", "SearchPoints" -> 100}]

it computes about several hours. The problem get worse if I consider higher order terms like

wG[t] = s0 Integrate[G[t1] HeavisideTheta[t - t1],{t, 0, t}] + s1 Integrate[G[t1] t1 HeavisideTheta[t - t1],{t1, 0, t}]

I suspect that this is because of ParametricNDSolve, since when I use NDSolve with fixed values of s, NMinimize works well.

So, my question is: is there a way of fast minimization of the distance with respect to s and s0 simultaneously?. Thank you very much in advance.

Some plots When I use 2 terms for wG[t], the two solutions are as follows:

while 3 terms give the following plot:

But to get this I need to run the code for over 30 times for different values of s.

Answer 1

A more or less clean approach, which relies on interpolating functions, Interpolation:

Clear[w, G, wG, distance]
w = NDSolveValue[{w'''[t] + 6 w[t] w'[t] - w'[t] == UnitStep[t], w[0] == w'[0] == w''[0] == 0}, w, {t, 0, 5}];
wG = ParametricNDSolveValue[{G'''[t] + 6 G[t] G'[t] - G'[t] == 0, G[0] == G'[0] == 0, G''[0] == s, wg'[t] == s0 G[t], wg[0] == 0}, {wg, G}, {t, 0, 5}, {s, s0}];
distance = Interpolation[Flatten[Table[{s, s0, NMaximize[{Abs[w[t] - wG[s, s0][[1]][t]], 0 <= t <= 5}, t][[1]]}, {s, 0, 3, .5}, {s0, 0, 5, .5}], 1]];
NMinimize[{distance[s, s0], 0 < s < 3, 0 < s0 < 5}, {s, s0}]
(* {0.388008, {s -> 2.14084, s0 -> 0.408041}} *)

It seems you are measuring distances with the $L^\infty$ norm; if you were to content yourself with an $L^1$ or $L^2$ norm, I guess one could speed up the computation a bit. If you want to consider more complex problems, then perhaps it will be necessary to do so, but I don't know whether such norms are acceptable for your problem.