5
$\begingroup$

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: <code>wsol[t]</code> and <code>wG[t]</code>

while 3 terms give the following plot:

enter image description here

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

$\endgroup$
2
$\begingroup$

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}} *)

enter image description here

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.

$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – Asatur Khurshudyan Jun 21 '18 at 16:18
  • $\begingroup$ @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. $\endgroup$ – AccidentalFourierTransform Jun 21 '18 at 16:20
  • $\begingroup$ 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? $\endgroup$ – Asatur Khurshudyan Jun 21 '18 at 16:29
  • 1
    $\begingroup$ @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. $\endgroup$ – AccidentalFourierTransform Jun 21 '18 at 16:51
  • 1
    $\begingroup$ @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. $\endgroup$ – AccidentalFourierTransform Jun 21 '18 at 23:01

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.