# Solving the non-linear heat equation numerically

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?

• 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
• 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: mathematica.stackexchange.com/q/73961/1871 . – xzczd Nov 30 '16 at 11:19
• 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
• 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