# Solving PDE with power series

I would like to solve the PDE

$$\partial_{x}f(x,y) + f(x,y)^2 = g(x,y)$$

with $$f(0,0)=0$$ and $$\partial_y f(0,0)=0$$ using a power series ansatz, i.e. I have an explicit expression for $$g(x,y)=\sin(x+y)\cos(x)\cos(2y)$$ and want to find the first coefficients of a power series $$f(x,y) = \sum_{n} c_n x^{n_1}y^{n_2}$$ solving this equation. I wonder if there is a way in mathematica to do this to get the first let's say $$\vert n\vert\le 11$$ coefficients.

• Use Series to expand the D[f[x,y],x]+f[x,y]^2-V[x,y]about {0, 0} and Solve for the coefficients. Apr 5, 2023 at 4:31
• @bbgodfrey I am very sorry, but since I am not very familiar with Mathematica, do you think you could provide a minimal working example? Apr 5, 2023 at 4:39
• @bbgodfrey In particular, how do I implement the initial conditions? Apr 5, 2023 at 4:45
• The expansion I suggested will contain f[0,0] and f[0,y]/.y->0. Set them to zero. First try the expansion with a small n. Apr 5, 2023 at 4:52
• @RolandF this cannot be solved by DSolve. Apr 5, 2023 at 11:29

With the additional initial conditions specified in a comment, the power series can be obtained by expanding the PDE as a series and solving the resulting algebraic equations. For instance,

n = 3;
CoefficientList[Normal@Series[
D[f[x, y], x] + f[x, y]^2 - Sin[x + y] Cos[x] Cos[2 y], {x, 0, n}, {y, 0, n + 1}]
/. {f[0, 0] -> 0, Derivative[0, _][f][0, 0] -> 0}, {x, y}];
Union@Flatten[%];
Solve[Thread[% == 0], Union@Variables@%] // Flatten;
Sum[Derivative[m1, m2][f][0, 0] x^m1 y^m2/{m1! m2!}, {m1, 1, n + 1},
{m2, 0, n + 1}] /. %

(* {x^2/2 - x^4/6 + x y - x^3 y/3 - x^4 y/4 - 5 x^2 y^2/4 -
x^3 y^2/3 + 5 x^4 y^2/12 - 13 x y^3/6 + 13 x^3 y^3/18 +
7 x^4 y^3/6 + 41 x^2 y^4/48 + 13 x^3 y^4/9 - 41 x^4 y^4/144} *)


The results rapidly becomes larger as n is increased.

• Problem 1: This is good in principle, but not so good in practical terms of how far it gives a decent approximation. I'll show a different approach and you can compare the discrepancies as (x,y) increase. Problem 2: I meant to upvote anyway, but apparently hit the wrong arrow and now cannot change it. Please edit so I can replace the errant downvote. Apr 5, 2023 at 15:23
• @DanielLichtblau I agree that a power series is not the optimal approach to solving this PDE, but perhaps the OP had a specific reason for requesting it. Also, I have edited my answer, as you requested. Apr 5, 2023 at 16:32
• Yes, maybe there was a reason. Thanks for editing, vote problem solved. (They should all be so simple.) Apr 5, 2023 at 16:46

Rather than use a series expansion, you might treat this as a parametric ODE and use ParametricNDSolve. Here of course y is the parameter, as it does not appear in the derivatives. I just solve from x=0 to x=2.

psol = ParametricNDSolveValue[{D[f[x], x] + f[x]^2 ==
Sin[x + y]*Cos[x]*Cos[2*y], f[0] == 0}, f[x], {x, 0, 2}, y];


Define an evaluation function.

pvalue[x0_, y0_] := psol[y0] /. x -> x0


I'll show this for a fairly modest range.

xmax = 1;
ymax = .7;
p1 = Plot3D[pvalue[x, y], {x, 0, xmax}, {y, 0., ymax}]


To compare to the series solution proposed by @bgodfrey we can do as follows.

pxy = x^2/2 - x^4/6 + x y - x^3 y/3 - x^4 y/4 - 5 x^2 y^2/4 -
x^3 y^2/3 + 5 x^4 y^2/12 - 13 x y^3/6 + 13 x^3 y^3/18 +
7 x^4 y^3/6 + 41 x^2 y^4/48 + 13 x^3 y^4/9 - 41 x^4 y^4/144;

p2 = Plot3D[pxy, {x, 0, xmax}, {y, 0., ymax}]


One can see these are close for a while but start to disagree, somewhat violently, in one corner (could use Show[{p1,p2}] to see it in more detail).

I should add that this method too will have its limitations-- the plot really does get wild beyond y=1.2 or so. I do not know if it's intrinsic to the problem or rather a limitation in the capabilities of ParametricNDSolve. Some experimenting also suggests it has no issue with y<0 until somewhere past y=-4.5.

• Problem 3: I munged the adding of images. How do I fix this? Apr 5, 2023 at 15:56
• Images now visible. Apr 5, 2023 at 16:26
• @bgodfrey What was the magic word? Apr 5, 2023 at 16:27
• I pressed <Enter> immediately after "past y=-4.5.". Apr 5, 2023 at 16:35
V[x_, y_] := Sin[x + y] Cos[x] Cos[2 y]
f[x_, y_] := Sum[c[n1, n2] x^n1 y^n2, {n1, 0, 11}, {n2, 0, 11}]
L = D[#, x] + #^2 &;
eqn = L[f[x, y]] - V[x, y] == 0;
sol = Solve[{f[0, 0] == 0, D[f[x, y], y] /. y -> 0 == 0}, Flatten@Table[c[n1, n2], {n1, 0, 11}, {n2, 0, 11}]];
Table[c[n1, n2] /. sol[[1]], {n1, 0, 2}, {n2, 0, 2}]
Expand[eqn /. sol[[1]]] // Simplify

• I reduced 11 to 1 and still this code just keeps on running. Is there a way to make it faster? Apr 5, 2023 at 4:02
• Your code contains errors. Also, may I suggest that, when you provide an answer, please include a sample result from running your code. Apr 5, 2023 at 4:45

If I have not made an error of reading the question, its a trival DSolve problem, not really a partial differential equation. We solve without that cos 2y term for a constant a. With cos y -term the equation is returned unsolved -our known complexity problem with trigonometry trees in the forests.

sol[x_, y_] =  f[x] /. DSolve[f '[x] ==  a  Sin[x + y] Cos[x] - f[x]^2, f,
x] [[1]] /. {a -> Cos[2 y]};

sol1[x_,y_]=sol[x,y]/.Solve[sol[0,0]==0,C[1]][[1]];

D[sol1[x, y], x] + sol1[x, y]^2 // Simplify;

Out=Cos[x] Cos[2 y] Sin[x + y]