1

I need to solve the following equation $\frac{\partial^2 \Phi}{\partial u^2}=\frac{\partial \Phi}{\partial u} \frac{\partial \Phi}{\partial \tau}$ with some initial condition $\Phi (u,0)$ to which the non-linear heat equation $\frac{\partial \Omega}{\partial \tau}=\frac{\partial^2 \mathrm{Log} \ \Omega}{\partial u^2}$ can be reduced by the substitution $\Omega=\frac{\partial \Phi}{\partial u}$.

Mathematica doesn't want to solve the original equation at all showing various singularities, but the reduced one is solved using NDSolve with very good accuracy (I checked it by direct substitution of the numerical solution in the equation).

But when I try to check if the original equation is also satisfied, Mathematica shows a terrible discrepancy. By various manipulations with the equation I discovered that the problems occur when there are $\frac{\partial^2 \Phi}{\partial u \ \partial \tau}$, $\frac{\partial^2 \Phi}{\partial \tau^2}$ or $\frac{\partial^3 \Phi}{\partial u^2 \partial \tau}$ terms in the expression I check, while there are no problems with $\frac{\partial \Phi}{\partial \tau}$, $\frac{\partial \Phi}{\partial u}$, $\frac{\partial^2 \Phi}{\partial u^2}$ and $\frac{\partial^3 \Phi}{\partial u^3}$ terms.

The code is the following:

Sol=NDSolve[{D[a[u, t], t] D[a[u, t], u] == D[a[u, t], {u, 2}], a[u, 0] == 1/2 Sqrt[1 + E^(2 u)]}, a, {t, 0, 2^35}, {u, -4.5, 20}, PrecisionGoal -> \[Infinity]]

Then I check numerically that the solution is OK:

Max[Abs[Table[N[((D[a[u, t], t] D[a[u, t], u])/D[a[u, t], {u, 2}] - 1) /. Sol /. u -> n/10 /. t -> 2^(k/10) - 1], {k, 0, 350}, {n, 20, 200}]]]

With the output:

0.000462987

So the precision is very good. But when I substitute it to the original equation and check the accuracy there:

Max[Abs[Table[N[(D[a[u, t], u, t]/D[Log[D[a[u, t], u]], {u, 2}] - 1) /. Sol /. u -> n/10 /. t -> 2^(k/10) - 1], {k, 0, 350}, {n, 20, 200}]]]

The output is:

644.342

So the accuracy is terrible.

The question is the following: does anybody understand why does NDSolve have problems with these terms, and how these problems could be fixed?

Andrey Feldman
  • 111
  • 3
  • 1
    Because substituting $\Omega=\frac{\partial \Phi}{\partial u}$ to $\frac{\partial \Omega}{\partial \tau}=\frac{\partial^2 \mathrm{Log} \ \Omega}{\partial u^2}$ won't produce $\frac{\partial^2 \Phi}{\partial u^2}=\frac{\partial \Phi}{\partial u} \frac{\partial \Phi}{\partial \tau}$?: With[{omega = Derivative[0, 1][phi][t, u]}, D[omega, t] == D[Log@omega, u, u]] – xzczd Nov 30 '16 at 11:04
  • @xzczd It does. Substitute it and integrate over $u$. An arbitrary function of $t$, which appears after the integration, is fixed by the appropriate choice of initial condition, and I set it to 0. – Andrey Feldman Nov 30 '16 at 11:07
  • 1
    Are you sure $\Phi$ is smooth enough in the whole domain? This is necessary condition for the equivalence of With[{omega = Derivative[0, 1][phi][t, u]}, D[omega, t] == D[Log@omega, u, u]] and 0==Integrate[Subtract@@With[{omega = Derivative[0, 1][phi][t, u]}, D[omega, t] == D[Log@omega, u, u]], u]. If you're sure, then the b.c. should be also added in NDSolve. Letting NDSolve add an artificial boundary condition is very dangerous, because no one actually knows what's added. Check this post: http://mathematica.stackexchange.com/q/73961/1871 . – xzczd Nov 30 '16 at 11:19
  • Even $\Omega$ looks very smooth and nice: Plot[Evaluate[Table[Log[D[a[u, 2^k - 1], u] /. Sol], {k, 0, 35}]], {u, 2, 20}, PlotRange -> All, ImageSize -> Full]. As for boundary conditions -- I specified them on the $\tau=0$ slice. Perhaps there are some better algorithms to solve PDEs than NDSolve? – Andrey Feldman Nov 30 '16 at 11:34
  • 1
    Don't you get a warning NDSolve::bcart that there aren't enough BCs? It expects one for a boundary like u == -4.5. (NDSolve calls the condition at t == 0 for the time integration an "initial condition" to distinguish between spatial and time BCs.) – Michael E2 Nov 30 '16 at 11:38
  • 1
    That's just an initial condition, you still need to add 2 boundary condition in $u$ direction. Or you're solving an initial value problem i.e. $-\infty<u<\infty$ i.e. there's no explicitly given boundary condition? – xzczd Nov 30 '16 at 11:40
  • @MichaelE2 You are right. I got this warning, but even with it, the reduced equation is solved with great accuracy as I pointed out, but the problems occured with checking its modified form. I had to restricted the entire real line to just an interval [-4.5,20] because there is no way to solve the equation on the infinite line. Do you know a way to circumvent this problem? – Andrey Feldman Nov 30 '16 at 11:43

0 Answers0