Numerical solution of KdV equation

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.

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.

• Thank you for your approach. I will check it now. Is there a way to consider also higher order terms of the form: wG[t] = Sum[s[n] Integrate[t1^n g[s][t1]] UnitStep[# - t1],{n, 0, N}] for a given N and minimize the distance with respect to s[n] simultaneously? Regarding the choice of norm. Actually, any norm can be used as soon as the two solutions get closer at least graphically. – Asatur Khurshudyan Jun 21 '18 at 16:18
• @AsaturKhurshudyan Yes, it is rather easy: just replace s0 Integrate[...] by s0 Integrate[...] + s1 Integrate[ t1 ...] + ... with as many terms as you wish. Note that the complexity grows fast with the number of parameters $s_i$, so it will take a lot of time. – AccidentalFourierTransform Jun 21 '18 at 16:20
• I replaced the constraint 0 < s0 < 1 by -10 < s0 < 10 (1 term) and it is still running. Maybe I should replace Integrate[...] by NIntegrate[...] or by Simpson's rule? – Asatur Khurshudyan Jun 21 '18 at 16:29
• @AsaturKhurshudyan Indeed, it takes more time, but replacing Integrate will not probably help (it is a rather fast operation, acting on an interpolating function, and it is essentially equivalent to NIntegrate). I'll check if memoization can speed things up a bit, and I'll let you know. – AccidentalFourierTransform Jun 21 '18 at 16:51
• @AsaturKhurshudyan You have a point, but I was about to leave. I'll check again later on, and if I can simplify the code further, I'll update the post. – AccidentalFourierTransform Jun 21 '18 at 23:01