# Accurate NDSolve for PDE with NIntegrate Functions

I have a function $$A(c)$$ which is defined by an integral. I wasn't able to solve it analytically, and therefore defined $$A(c)$$ using NIntegrate:

A[c_?NumericQ] := 2 NIntegrate[Sqrt[(3 + 5 y^4 + c (1 + 4 y^2 - y^4) - Sqrt[(-1 + c) (-9 + 42 y^2 - 17 y^4 + c (-7 - 10 y^2 + y^4))] + y^2 (-12 + Sqrt[(-1 + c) (-9 + 42 y^2 - 17 y^4 + c (-7 - 10 y^2 + y^4))]))/(1 + c + (-5 + c) y^2 + 2 y^4)]/Sqrt[2],
{y, -Sqrt[((3 - Sqrt[9 - 7 c^2 - 2 c] + c)/(2 + 2 c))], Sqrt[(3 - Sqrt[9 - 7 c^2 - 2 c] + c)/(2 + 2 c)]}]


Using $$A(c)$$ I define a function $$rsqr(x,y)$$:

rsqr[x_, y_] := 1/(2 π)A[1 - 4/(1/(x + 1)^2 + 1/(y + 1)^2 + 1/(x - 1)^2 + 1/(y - 1)^2)]


I want to solve the following first order PDE for the function $$\alpha(x,y)$$: $$\partial_{x}rsqr\partial_{y}\alpha-\partial_{y}rsqr\partial_{x}\alpha=1$$ The region is $$[0,1]\times[0,1]$$ rectangle, and the boundry condition is $$\alpha(x,0)=0$$ (I know it is insufficient boundry conditions, but doesn't seem to be a problem). First, I have to calculate the derivatives of $$rsqr(x,y)$$. Following another question I have asked, I calculate the derivatives numerically using finite-differences:

h = 10^-5;
drsqrdx[x_, y_] := Piecewise[{{(rsqr[x + h, y] - rsqr[x - h, y])/(2 h), h < x < 1 - h},
{(rsqr[x, y] - rsqr[x - h, y])/h, x >= 1 - h},
{(rsqr[x + h, y] - rsqr[x, y])/h, x <= h}}]
drsqrdy[x_, y_] := Piecewise[{{(rsqr[x, y + h] - rsqr[x, y - h])/(2 h), h < y < 1 - h},
{(rsqr[x, y] - rsqr[x, y - h])/h, y >= 1 - h},
{(rsqr[x, y + h] - rsqr[x, y])/h, y <= h}}]


Now I try to solve the PDE:

NDSolve[{drsqrdx[x, y] \!$$\*SubscriptBox[\(∂$$, $$y$$]$$α[x, y]$$\) - drsqrdy[x, y] \!$$\*SubscriptBox[\(∂$$, $$x$$]$$α[x, y]$$\) == 1, α[x, 0] == 0}
,α, {x, 0, 1}, {y, 0, 1}]


But for some reason it gives a warning: Power::infy: Infinite expression 1/0.^2 encountered.. Therefore I try to solve in $$[0,1-\epsilon]\times[0,1-\epsilon]$$ rectangle for small $$\epsilon$$ in hope the interpolating function would approximate it in the whole rectangle:

ϵ = 10^-6;
αsol = NDSolve[{drsqrdx[x, y] \!$$\*SubscriptBox[\(∂$$, $$y$$]$$α[x, y]$$\) - drsqrdy[x, y] \!$$\*SubscriptBox[\(∂$$, $$x$$]$$α[x, y]$$\) == 1, α[x, 0] == 0}
,α, {x, 0, 1-ϵ}, {y, 0, 1-ϵ}]


It takes a long time and produces warnings, but gives a solution. The problem is that the solution does not satisfy the equation. By plotting $$\partial_{x}rsqr\partial_{y}\alpha-\partial_{y}rsqr\partial_{x}\alpha-1$$ we can see a significant deviation from $$0$$ in some places, especially near the origin:

dαdx[x_, y_] := Piecewise[{{(α[x + h, y] - α[x - h, y])/(2 h), h < x < 1 - h},
{(α[x, y] - α[x - h, y])/h, x >= 1 - h},
{(α[x + h, y] - α[x, y])/h, x <= h}}] /. αsol[[1]]
dαdy[x_, y_] := Piecewise[{{(α[x, y + h] - α[x, y - h])/(2 h), h < y < 1 - h},
{(α[x, y] -α[x, y - h])/h, y >= 1 - h},
{(α[x, y + h] - α[x, y])/h, y <= h}}] /. αsol[[1]]
Plot3D[{drsqrdx[x, y] dαdy[x, y] - drsqrdy[x, y] dαdx[x, y] - 1, 0},
{x, 0, 1}, {y, 0, 1}, AxesLabel -> Automatic, PlotRange -> {-0.5, 0.5}]


I tried to change many options, such as adding MaxStepSize -> 0.01 on NDSolve, increasing AccuracyGoal, PrecisionGoal or WorkingPrecision for NDSolve and NIntegrate, and trying many Methods of NDSolve and NIntegrate. While it does give better results in some cases, and does not produce warnings (except insufficient boundry conditions), the results still never satisfy the PDE.

Does anyone know how to fix this problem and find a solution to the PDE? Thank you!

Based on the Knowledge of A[c]  here we try to create a FEM-solution using ElementMeshInterpolation.

mesh definition:

Needs["NDSolveFEM"]
mesh = ToElementMesh[Rectangle[{0, 0}, {1 , 1 }]   ,"MeshElementType" ->TriangleElement  , "MeshOrder" -> 1 , "MaxCellMeasure" -> .001  ]
mesh["Wireframe"]


elementmeshinterpolation:

xy = mesh["Coordinates"];
z = Map[1/(2 Pi) A[arg[ #[[1]], #[[2]] ] ]  &, xy] ;
rsqr = ElementMeshInterpolation[mesh, z ]; (* rsqr[x,y]*)


ElementMeshInterpolation allows extrapolation, that's why we don't need piecewise-functionality

h = 10^-3;
drsqrdx = ElementMeshInterpolation[mesh,Map[(rsqr[#[[1]] + h, #[[2]]] - rsqr[#[[1]] - h, #[[2]]])/(2 h) &, xy]];
drsqrdy = ElementMeshInterpolation[mesh,Map[(rsqr[#[[1]], #[[2]] + h] -rsqr[#[[1]], #[[2]] - h])/(2 h) &,  xy]];


FEM solution:

pde = drsqrdx[x, y] Derivative[0, 1][\[Alpha]][x, y] -drsqrdy[x, y] Derivative[1, 0][\[Alpha]][x, y] == 1

alfa = NDSolveValue[{pde, \[Alpha][x, 0] == 0}, \[Alpha],Element[{x, y}, Rectangle[{0, 0}, {1 , 1 }]  ], AccuracyGoal -> 15,PrecisionGoal -> 8]
Plot3D[alfa[x, y], Element[{x, y}, Rectangle[{0, 0}, {1 , 1 }]  ],AxesLabel -> {x, y, \[Alpha]}]


"approximation error":

res = (pde /. Equal -> Subtract) /. \[Alpha] -> alfa;
Plot3D[res, Element[{x, y}, mesh], PlotRange -> All,AxesLabel -> {"x", "y", "res"} ]


• Thank you again @Ulrich. I have checked this method and it seems to be very helpful. The error you show is still pretty large. By using the options "MeshOrder" -> 2 and "MaxCellMeasure" -> .00001 it got better, but it has points where it is still large (up to about 0.5). Is there a way to make the error small, lets say up to $10^-3$? Taking smaller MaxCellMeasure doesn't seem to solve it, and takes a very long time to compute.
– Roy
Commented Dec 28, 2022 at 13:33
• @Roy FEM minimizesbweighted residuals of the pde. That's why substituting the solution into pde gives no error measure I think Commented Dec 28, 2022 at 14:16
• I don't think I understand what you are tring to say here, excuse me. Is there a way to minimize it "better", such that the error would be actually small?
– Roy
Commented Dec 28, 2022 at 18:48
• My point is: FEM solution fullfills the pde only for the limiting case of vanishing elementsize. But provides a good solution for small elements Commented Dec 28, 2022 at 20:19
• @Roy In link is an example of a pde problem which is solved with FEM and shows an accurate approximation compared to the known exact solution. If you plot the residuals Plot3D[pde - f /. u -> ufun, Element[{x, y}, r]] you'll see large residuals! Commented Jan 5, 2023 at 18:21

The main problem seems to be: You have a smooth approximation A[c], but your pde needs smooth derivative A'[c]! Perhaps you know this derivativefunction???

Because I don't have further knowledge about A'[c] I'll use an approximation A[c]~4 (1 - Sqrt[1 - c])

Plot[{A[c] , 4 (1 - Sqrt[1 - c]) }, {c, 0, 1},PlotStyle ->{Automatic, {Dashed, Red}}]


which fits quite well. Using

As = Apply[Function, {c, D[(4 (1 - Sqrt[1 - c]) ), c]}]


and

arg = Function[{x, y}, 1 - 4/(1/(x + 1)^2 + 1/(y + 1)^2 + 1/(x - 1)^2 + 1/(y - 1)^2)]    (* function argument inside rsqr[x,y] ...*)


we get the solution

pde = 1/(2 Pi)As[arg[x, y]] Derivative[1, 0 ][arg][x,y] Derivative[0, 1 ][\[Alpha]][x, y] -
1/(2 Pi) As[arg[x, y]] Derivative[0, 1  ][arg][x,y] Derivative[1, 0  ][\[Alpha]][x, y] == 1
alfa = NDSolveValue[{pde, \[Alpha][x, 0] == 0}, \[Alpha],Element[{x, y}, Rectangle[]] ]
Plot3D[alfa[x, y], Element[{x, y}, Rectangle[]]]


Hope it helps!

• Dear @Ulrich, thank you for your answer! 1. Unfortunately I do not know much about the derivative of A(c). How did you know that 4(1-Sqrt[1-c]) would be such a good approximation? While this approximation is interesting, it is not good enough for what I'm trying to do, so I'll have to stick with A(c) itself. 2. I don't think the problem is the smoothness of the derivatives. The solution you show here also have the same problem I showed in the original post: by plotting you can verify that it does not satisfy the PDE. The solution I got using the full A(c) looks almost the same.
– Roy
Commented Dec 26, 2022 at 22:20
• @Roy #1: You might improve the approximation by fitting  A[c]~(4 (1 - Sqrt[1 - c]) + c1 c Sqrt[1 - c]^3. #2: Numerical differentiation is usually a problem. Be aware that the solution you got is a FEM-solution , which only gives good approximation at grid-points of the mesh used by NDSolve Commented Dec 27, 2022 at 10:42
• Thank you @Ulrich. It is a FEM-solution, but isn't the interpolating function gives a good approximation also between the grid points? Taking a denser grid doesn't seem to give a solution that satisfies the PDE near the origin. I'm aware that numerical differentiation is problematic, but is there a way to find a solution that satisfies the PDE if I use the approximation for A(c)? The solution you got also does not satisfy it.
– Roy
Commented Dec 27, 2022 at 11:45
• @Roy At the origin solution seems to jump for x=0,y=0.0001 . I' ll append a solution idea based on your approach using FEM& ElementMeshInterpolation Commented Dec 27, 2022 at 12:04