Numerically solve PDEs with constraints and without boundary solution

I have a PDE like

D[h[x1, x2], x1]*a[x1,x2]+D[h[x1,x2], x2]*b[x1,x2] + c[x1,x2] == h[x1,x2]


where a,b,c are known functions of x1 and x2, and h are the function to be solved. x1 and x2 are both in [-2, 2]. For some selected a,b,c, DSolveValue can give me perfect analytical solutions, like (a=0.1*x1, b=x2-x1^2, c=-x1^2)

eq1 = -x1^2 + D[h1[x1, x2], x1]*(0.1)*x1 +
D[h1[x1, x2], x2]*(x2 - x1^2) == h1[x1, x2];
s1 = DSolveValue[eq1, h1[x1, x2], {x1, x2}]
(*x1^2. (-1.25 + x1^8. C[1][(0.25 (-5. x1^2 + 4. x2))/x1^10])*)


For this analytical solution, I then use the gradient constraints to figure out how C[1] should look like and finally get h1=-1.25*x1^2, which is the solution I want. But for other a,b,c, DSolveValue reports that "Inverse functions are being used by Solve, so solutions may not be found". For example (a=x2, b=Sin[x1]+x2, c=0):

eq2 = D[h2[x1, x2], x1]*x2 + D[h2[x1, x2], x2]*(Sin[x1] + x2) ==
h2[x1, x2];
s2 = DSolveValue[eq2, h2[x1, x2], {x1, x2}]


The output just repeats my cmd and doesn't have a solution. I'm aware that this is because h2 doesn't have an analytical solution. How to get a numerical solution under this circumstance?

• First we need to restrict x1, x2 as well. Jan 21, 2021 at 13:50
• Your PDE does not define the derivatives at x1==x2==0 Jan 21, 2021 at 13:53

Linear PDEs typically are solved by the method of characteristics. For the PDE in the question, the ODEs of the characteristics are

{x2'[s] == Sin[x1[s]] + x2[s], x1'[s] == x2[s], h'[s] == h[s]}


Attempting to solve these ODEs with DSolve is unsuccessful, returning the Solve::ifun error message. That this is the same error cited in the question is not surprising, because DSolve itself undoubtedly uses the method of characteristics when attempting to solve a linear or quasilinear PDE. (See questions 238346 and 238808 for examples in which the characteristics, and therefore the PDE, can be solved symbolically.) However, NDSolve can evaluate the characteristics with little difficulty. Focus first on the ODEs for (x1, x2}.

sp = ParametricNDSolveValue[{x2'[s] == Sin[x1[s]] + x2[s], x1'[s] == x2[s],
x1[0] == x10, x2[0] == x20}, {x1[s], x2[s]}, {s, -10, 10}, {x10, x20}];
Table[sp[n, -2], {n, -4.03, 2, .1}];
pt = ParametricPlot[%, {s, -10, 10}, PlotRange -> {{-2, 2}, {-2, 2}},
ImageSize -> Large, FrameLabel -> {x1, x2}, LabelStyle -> {15, Bold, Black},
PlotStyle -> RGBColor[0.368417, 0.506779, 0.709798]];


Before displaying the result, let us superimpose (in black) the two stream lines that pass through the origin.

NDSolveValue[{x2'[s] == Sin[x1[s]] + x2[s], x1'[s] == x2[s], x1[10^-5] == 10^-5,
x2[10^-5] == 10^-5 1/2 (1 - Sqrt[5])}, {x1[s], x2[s]}, {s, -50, 50}];
pt0 = ParametricPlot[%, {s, -50, 50}, PlotStyle -> Black,
PlotRange -> {{-2, 2}, {-2, 2}}];
Show[pt, pt /. Line[z_] -> Line[-z], pt0, pt0 /. Line[z_] -> Line[-z]]


(The substitution, Line[z_] -> Line[-z] reflects the solution in the origin, reducing computations by a factor of two.) Now, the point of evaluating the characteristics is that h2[x1, x2] satisfies

DSolveValue[h'[s] == h[s], h[s], s]
(* E^s C[1] *)


along each stream line, with C[1] a constant that can vary from stream line to stream line. Consequently, if the value of h2[x1, x2] is known anywhere on each stream line, its value is known everywhere. The question above specifies only that the gradient of h2[x1, x2] vanishes at the origin. From the PDE, it immediately follows that h2[0, 0] also vanishes, and the value of h2[x1, x2] must be zero everywhere along the black lines. Moreover, h[x1, x2] also must vanish on characteristics passing infinitesimally close to the origin. However, values of h2[x1, x2] on other stream lines are unspecified.

For completeness, here is a plot of the stream lines over a larger domain.

rl = Table[{Line[z_] :> Line[{2 nn Pi + #[[1]], #[[2]]} & /@ z]} /.
nn -> n, {n, -2, 2}];
Table[sp[n, -4], {n, -6 Pi, 3 Pi, Pi/5}];
ParametricPlot[%, {s, -10, 10}, PlotRange -> {{-10, 10}, {-4, 4}},
ImageSize -> Large, FrameLabel -> {x1, x2}, LabelStyle -> {15, Bold, Black},
PlotStyle -> RGBColor[0.368417, 0.506779, 0.709798]];
Show[%, % /. Line[z_] -> Line[-z], Replace[{pt0, pt0 /. Line[z_] -> Line[-z]},
rl, All]]


As would be expected, the structure of the stream lines is periodic in x1 with period 2 Pi. Finally, it should be noted that StreamPlot can produce similar plots, although with less detail.

Response to OP's comment

Le ZHENG asked in a comment below about the initial conditions in

    NDSolveValue[{x2'[s] == Sin[x1[s]] + x2[s], x1'[s] == x2[s], x1[10^-5] == 10^-5,
x2[10^-5] == 10^-5 1/2 (1 - Sqrt[5])}, {x1[s], x2[s]}, {s, -50, 50}];


used above. They were obtained by combining the first two characteristic equations at the beginning of the answer to obtain

x2'[x1] == (Sin[x1] + x2[x1])/x2[x1]


For stream lines passing through the origin, x2'[x1] is equal to x2[x1]/x1. With these expressions set to a for convenience, the preceding equation becomes a == 1/a + 1, which has the solutions

Solve[a == 1/a + 1, a] // Flatten
(* {a -> 1/2 (1 - Sqrt[5]), a -> 1/2 (1 + Sqrt[5])} *)


Hence, very near the origin x2 = x1 1/2 (1 - Sqrt[5]) for one stream line and x2 = x1 1/2 (1 + Sqrt[5]) for the other. I chose "very near" to be x1 = 10^-5. And, because the set of characteristic equations is autonomous (not explicitly dependent on s), I could choose s corresponding to x1 = 10^-5 to be any value whatsoever. I chose it equal to x1.

By the way, no stream lines pass through {2 PI n, 0}. with n an odd integer. Instead, they spiral around those points with exponentially decreasing radii.

Sample Complete Numerical Solution

As a sample solution, choose h[x1,x2] = (x1^2 +x2^2)/4 along the x and y axes, which is sufficient to completely determine h and satisfies the requirements that h and its gradients vanish at the origin. It is convenient for this computation to eliminate s from the characteristic ODEs in favor of x1 or x2.

sp1 = ParametricNDSolveValue[{x2'[x1] == (Sin[x1] + x2[x1])/x2[x1],
h'[x1] == h[x1]/x2[x1], x2[0] == x0, h[0] == x0^2/4,
WhenEvent[x2[x1] > 3.5, "StopIntegration"]}, {x1, x2[x1], h[x1]}, {x1, -2, 2}, x0];
sp2 = ParametricNDSolveValue[{x1'[x2] == x2/(Sin[x1[x2]] + x2),
h'[x2] == h[x2]/(Sin[x1[x2]] + x2), x1[0] == x0, h[0] == x0^2/4,
WhenEvent[x1[x2] > 2, "StopIntegration"]}, {x1[x2], x2, h[x2]}, {x2, -2, 2}, x0];

plt1 = Show@Table[tem = sp1[n]; ParametricPlot3D[tem,
Evaluate[Join[{x1}, Flatten[Last[tem] /. x1 -> "Domain"]]],
PlotRange -> {{-2, 2}, {-2, 2}, {0, 1.5}}, AxesLabel -> {x1, x2, h},
ImageSize -> Large, LabelStyle -> {15, Bold, Black}], {n, .01, 3.5, .05}];
plt2 = Show@Table[tem = sp2[n]; ParametricPlot3D[tem,
Evaluate[Join[{x2}, Flatten[Last[tem] /. x2 -> "Domain"]]],
PlotRange -> {{-2, 2}, {-2, 2}, {0, 1.5}}, AxesLabel -> {x1, x2, h},
ImageSize -> Large, LabelStyle -> {15, Bold, Black}], {n, .01, 2, .05}];
Show[plt1, plt1 /. Line[z_] :> Line[Map[{-#[[1]], -#[[2]], #[[3]]} &, z]],
plt2, plt2 /. Line[z_] :> Line[Map[{-#[[1]], -#[[2]], #[[3]]} &, z]]]


As expected, h = 0 on the stream lines passing through the origin and those infinitesimally close. The solution also can be displayed with Plot3D.

Flatten[Cases[%, Line[z_] -> z, Infinity], 1];
ListPlot3D[%, PlotRange -> {{-2, 2}, {-2, 2}, {0, 1.5}}, AxesLabel -> {x1, x2, h},
ImageSize -> Large, LabelStyle -> {15, Bold, Black}]


• The answer is fantasy. Let me try to digest it. Thank you very much. Jan 26, 2021 at 0:12
• hi bbgodfrey, where does x1[10^-5] == 10^-5, x2[10^-5] == 10^-5 1/2 (1 - Sqrt[5]) come from in the code? Jan 30, 2021 at 10:27
• @LeZHENG I have addressed you last comment at the end of my answer above. By the way, I suggest that you 1) take the introductory tour now, 2) when you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge, 3) remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign, and 4) give help too, by answering questions in your areas of expertise. Jan 30, 2021 at 21:48

It also is possible to obtain an approximate symbolic solution, valid for small x1. Expand the PDE in x1 to obtain

DSolve[(x1 + x2)*D[h[x1, x2] + x2*D[h[x1, x2] == h[x1, x2], h[x1, x2], {x1, x2}]


Although DSolve cannot solve this PDE directly, it can solve the characteristic ODEs.

{x2'[s] == x1[s] + x2[s], x1'[s] == x2[s], h'[s] == h[s]}


Focus first on the ODEs for {x1, x2}

DSolve[{x1'[s] == x2[s], x2'[s] == x1[s] + x2[s]}, {x1[s], x2[s]}, s] // Flatten
sol = ((Simplify[% /. C[n_] -> C[n] Exp[-s/2]] // ExpToTrig // Simplify)
/. C[n_] -> C[n] Exp[s/2]) // Simplify
(* {x1[s] -> 1/5 E^(s/2) (5 C[1] Cosh[(Sqrt[5] s)/2] -
Sqrt[5] (C[1] - 2 C[2]) Sinh[(Sqrt[5] s)/2]),
x2[s] -> 1/5 E^(s/2) (5 C[2] Cosh[(Sqrt[5] s)/2] +
Sqrt[5] (2 C[1] + C[2]) Sinh[(Sqrt[5] s)/2])} *)


which can be plotted as

Show[ParametricPlot[Table[{x1[s], x2[s]} /. sol /. {C[1] -> n, C[2] -> 0},
{n, -2,   2, .1}], {s, -3, 3}, PlotRange -> {{-2, 2}, {-2, 2}},
ImageSize -> Large, LabelStyle -> {15, Bold, Black}],
ParametricPlot[Table[{x1[s], x2[s]} /. solt /. {C[1] -> 0, C[2] -> n}, {n, -3, 3, .1}],
{s, -3, 3}]]


Compare this to the first plot in my strictly numerical answer above. Clearly, they are about the same for small x1. Next, obtain a symbolic solution for h. Consider the quadratic expression, (x1[s] + x2[s]/2)^2 - 5/4 x2[s]^2), which reduces to

Simplify[((x1[s] + x2[s]/2)^2 - 5/4 x2[s]^2) /. sol /. C[2] -> 0]
(* E^s C[1]^2 *)


and

DSolve[h'[s] == h[s], h[s], s] /. C[1] -> C[3]
(* {{h[s] -> E^s C[3]}} *)


E^s can be eliminated between these two results to yield

h[x1, x2] == Simplify[(x1 + x2/2)^2 - 5/4 x2^2] C[3]/C[1]^2
(* h[x1, x2] == ((x1^2 + x1 x2 - x2^2) C[3])/C[1]^2 *)


Similarly, assuming C[2] -> 0 leads to

(* h[x1, x2] == ((-x1^2 - x1 x2 + x2^2) C[3])/C[2]^2 *)


Choose as boundary conditions h[x1,x2] = (x1^2 +x2^2)/4 along the x1 and X2 axes, as done in the numerical solution above. Then, C[3]/C[2]^2 reduces to 1/4 on the x1 axis, and C[3]/C[1]^2 reduces to 1/4 on the x2 axis. Combining these expressions gives for h

h[x1, x2] = Abs[(x1^2 + x1 x2 - x2^2)/4];


which can be plotted as

Plot3D[Abs[(x1^2 + x1 x2 - x2^2)/4], {x1, -2, 2}, {x2, -2, 2},
PlotRange -> {0, 1.5}, AxesLabel -> {x1, x2, h}, ImageSize -> Large,
LabelStyle -> {15, Bold, Black}, PlotPoints -> 50]


which agrees well for small x1 with the final plot in the numerical answer above.