# how to fix this PDE numerical solution (NDSolve)

I have the following code:

ClearAll[x, t, u];
a = 10;
sol = NDSolve[{D[u[x, t], x] == D[u[x, t], t], u[x, 0] == x^2,
u[0, t] == t}, u, {x, -a, a}, {t, 0, 10}, MaxStepSize -> 0.01]
Plot3D[{u[x, t] /. sol}, {x, -a, a}, {t, 0, 10}]


This results in the following graph :

As you can see, the graph is not smooth as it should be. I've tried to fix it by adding MaxStepSize -> 0.01. This changes the graph somewhat, but doesn't really fix it as you can see.

How do I get a smooth graph? what am I doing wrong?

## 2 Answers

First, I'd like to point out the problem doesn't have a classical solution, which should be $C^1$ continuous and satisfy all the i.c. and b.c. in the whole interested domain. This can be easily proven by finding the general solution of the equation:

{eq, ic, bc} = {D[u[x, t], x] == D[u[x, t], t], u[x, 0] == x^2, u[0, t] == t}
generalsol = First@DSolve[eq, u, {x, t}]
(* {u -> Function[{x, t}, C[1][t + x]]} *)


and substituting it to i.c. and b.c.:

{ic, bc} /. generalsol
(* {C[1][x] == x^2, C[1][t] == t} *)


Clearly the only function relationship satisfying ic is Function[a, a^2], while the one satisfying bc is Function[a, a], so they can't be satisfied simultaneously, as long as we require the solution to be smooth.

Still, if we lower our requirement a bit i.e. allow solution that doesn't satisfy the equation and corresponding i.c. and b.c. in some region, we can find weak solution. (Well, to be honest I'm not sure if it's the correct terminology. I've asked it in Math.SE.)

One weak solution can be found by using Laplace transform on t (I'll use pdeSolveWithLaplaceTransform for the task):

(* Definition of pdeSolveWithLaplaceTransform isn't included in this post,
please find it in the post above. *)
asol[x_, t_] = pdeSolveWithLaplaceTransform[{eq, bc}, ic, u[x, t], t, x]
(* (t + x)^2 + (t - t^2 + x - 2 t x - x^2) HeavisideTheta[t + x] *)


It's the same as given by DSolve, as shown in zhk's answer.

Notice this solution does not satisfy the original equation at $t+x=0$ (HeavisideTheta isn't defined at 0), ic when $x>0$ and bc when $t<0$:

Plot[Subtract @@ ic /. u -> asol // Evaluate, {x, -10, 10}]


This isn't the only possible weak solution, if we treat x as time i.e.:

asol2[x_, t_] = pdeSolveWithLaplaceTransform[{eq, ic}, bc, u[x, t], x, t]
(* t + x - (t + x) HeavisideTheta[t + x] + (t + x)^2 HeavisideTheta[t + x] *)


This solution does not satisfy the original equation at $t+x=0$, ic when $x<0$ and bc when $t>0$.

As you see, when weak solution comes on stage, the problem doesn't have a unique solution. If neither of the 2 solutions above is the one you need, there must be a missing auxiliary conditions. As to auxiliary condition, here's an example:

Position of discontinuous coefficient influences the solution of PDE

Finally, I believe the reason why NDSolve has trouble in finding the numeric solution of the problem is, currently NDSolve doesn't have a strong enough support for weak solution problem. Here's just another few troublesome cases:

1D Euler equations (fluid dynamics) with NDSolve

Circumvent NDSolve::bdord: error for 1-D Euler Equations

Solving wave equation with singular initial conditions

• To the downvoter, I am interested in what was missing from my answer, would you please elaborate? I'm not trying to complain here, I'm just curious about what I can do to improve my answer. – xzczd May 8 '17 at 7:52
• Your statement about a weak solution is unbased. Up to en.wikipedia.org/wiki/Weak_solution , a weak solution must satisfy the condition after " then an integrable function u would be said to be a weak solution if" in the cited Wiki article. This condition is not verified by you. Also there are problems with IC. @xzczd: – user64494 May 8 '17 at 8:17
• @user64494 Then I think LaplaceTransform, which is essentially a integration, plays exactly the role of that "condition". As to the i.c., the definition in the wiki page doesn't have any requirement for it, and I think the "problem" with i.c. does not contradict the definition: if the weak solution may not satisfy the equation in some "region", then why can't it be on the boundary? – xzczd May 8 '17 at 8:37
• I prefer solid arguments supplied with codes and references over unbased words. @xzczd – user64494 May 8 '17 at 11:01
• @user64494 I've rephrased my answer to make it more prudent. – xzczd May 9 '17 at 16:02

Why not use DSolve if it is solving the PDE?

ClearAll[x, t, u];
a = 10;
sol = DSolve[{D[u[x, t], x] == D[u[x, t], t], u[x, 0] == x^2, u[0, t] == t}, u, {x, -a, a},
{t, 0, 10}]

Plot3D[{u[x, t] /. sol}, {x, -a, a}, {t, 0, 10}]


• The exact solution does not seem to be true. Look at the output of FullSimplify[{D[u[x, t], x] == D[u[x, t], t], u[x, 0] == x^2, u[0, t] == t} /. %, Assumptions -> t > 0] @zhk – user64494 May 7 '17 at 15:47
• See the executed code in pdf file through dropbox.com/s/pyag64abroxmdvo/Bug%20in%20DSolve.pdf?dl=0 and in nb file through dropbox.com/s/5sjhkkizstlqq27/Bug%20in%20DSolve.nb?dl=0 @zhk – user64494 May 7 '17 at 16:06
• But :sol = DSolve[{D[u[x, t], x] == D[u[x, t], t], u[x, 0] == x^2, u[0, t] == t}, u, {x, -10, 10}, {t, 0, 10}]; FullSimplify[{D[u[x, t], x] == D[u[x, t], t], u[x, 0] == x^2, u[0, t] == t} /. sol, Assumptions -> {t > 0, x < 0}] is {{True, True, True}} – Mariusz Iwaniuk May 7 '17 at 19:03