# Partial differential equation, Finite difference Method

For my research project I am trying to solve following partial differential equation in Mathematica

(V[r]*ψ[u, v])/4 + Derivative[1, 1][ψ][u, v] == 0


with boundary conditions

ψ[u,v0]==0; ψ[u0,v]==Exp[-((v - vc)^2/(2 σ^2))]; vc = 10; σ = 3


Here (u,v) are two independent variables and ψ is a dependent variable. Above two-dimensional equation can be integrated numerically, using the finite difference method (suggested in some paper) as

Subscript[ψ, N] = Subscript[ψ, E] - Subscript[ψ, S] + Subscript[ψ, W] -
1/32 V δu δv (-Subscript[u, E] - Subscript[u, N] +
Subscript[v, N] + Subscript[v, W]) (Subscript[ψ, E] + Subscript[ψ, W])


Here N,W,E and S forms a rectangular grid means If you make a grid in u (y-axis) and v (x axis) then N represents (u+δu, v+δv), E represents (u, v+δv), W represents (u+δu,v) and S represents (u, v) points on the grid. δu and δv are grid spacing. One can take δu=δv=0.1 (or small). r -> { rp to Infinity}, v -> { - Infinity to 40}, u -> { - 40 to Infinity}. u0 and v0 are the starting points on the grid. You can take them anywhere provided they are consistent with above mentioned interval.

V[r]= ((-(1/(50 r^3)) + 2 r + 1/(100 r^2 rp) + rp/r^2 + rp^3/r^2) (1 + 1/(
100 r^2) + r^2 - 1/(100 r rp) - rp/r - rp^3/r))/r


is a function of "r" rather than (u,v). r is related to u and v as,

u=t-f[r], v=t+f[r]


where

f[r]= 1/(2 KP) Log[r - rp] - 1/(2 KN) Log[r - rN] -
A Log[r^2 + r P + q] + (2 (B + A P))/
Sqrt[4 q - P^2] (ArcTan[(2 r + P)/Sqrt[4 q - P^2]] - Pi/2) - 20

A = (R^2 (rp + rN) (rp^2 + rN^2 + 2 rp rN + R^2))/(
2 (3 rp^2 + rN^2 + 2 rp rN + R^2) (3 rN^2 + rp^2 + 2 rp rN + R^2))

B = (R^2 (rp^2 + rN^2 + R^2) (rp^2 + rN^2 + rp rN + R^2))/((3 rp^2 +
rN^2 + 2 rp rN + R^2) (3 rN^2 + rp^2 + 2 rp rN + R^2))

KP = ((rp - rN) (3 rp^2 + rN^2 + 2 rp rN + R^2))/(2 R^2 rp^2)

KN = ((rp - rN) (3 rN^2 + rp^2 + 2 rp rN + R^2))/(2 R^2 rN^2)

rN = -(rp/3) - (2^(1/3) (3 rp^2 + 2 rp^4))/(
3 rp ((27 rp^2)/100 + 9 rp^4 + 7 rp^6 + Sqrt[
4 (3 rp^2 + 2 rp^4)^3 + ((27 rp^2)/100 + 9 rp^4 + 7 rp^6)^2])^(
1/3)) + ((27 rp^2)/100 + 9 rp^4 + 7 rp^6 + Sqrt[
4 (3 rp^2 + 2 rp^4)^3 + ((27 rp^2)/100 + 9 rp^4 + 7 rp^6)^2])^(
1/3)/(3 2^(1/3) rp)

q=rp^2 + rN^2 + rp rN + R^2

P=rp + rN

rp = 4/10

R=1.


So now the idea is for a particular value of u and v, f[r]=(v-u)/2 is fixed. Then we can use f[r] equation to find r and from this r we can find V[r]. This will give V for those particular values of u and v. With this one can solve the partial differential equation and find ψ[u,v].

I need a plot of ψ[400, v] as a function of v. The plot of ψ[400, v] will be oscillating and decaying. An image of the result is shown below. I would really appreciate if someone can help in this regard.

• Have you actually integrate the equation numerically, as you say is possible early in your question? Please show the calculations you actually have carried out. – bbgodfrey Nov 27 '15 at 6:03
• @bbgodfrey...Yes, I have tried to find [Psi][u,v] using Subscript[[Psi], N] = Subscript[[Psi], E] + Subscript[[Psi], W] - Subscript[[Psi], S] - ([Delta]u [Delta]v (Subscript[[Psi], E] + Subscript[[Psi], W]))/8 V (( Subscript[v, N] + Subscript[v, W] - Subscript[u, N] - Subscript[u, E])/4) , but could not reproduce the results. – mat Nov 27 '15 at 6:08
• Boundary conditions need to be specified in order to solve the first equation numerically. Also, your second equation has no obvious connection to the rest of your question, so it is not surprising that you could not compute ψ from it. – bbgodfrey Nov 27 '15 at 6:18
• R remains undefined, and the question now contains two inconsistent definitions for rN. Yu do not make it easy for people to help you. – bbgodfrey Nov 27 '15 at 14:26
• I managed to solve the equation, as soon as this question is reopened, I'll post an answer. But with my arbitrarily chosen u0 and v0, the image I obtained is quite different from yours. – xzczd Nov 28 '15 at 8:57

As mentioned in the comment above, I didn't reproduce the image you gave, maybe my arbitrarily chosen u0 and v0 are not proper, maybe I've made a mistake somewhere, maybe something is wrong with the equation you provided. Anyway I'll show my 2 solutions below, one with NDSolve, the other with FDM using the difference scheme you provided, the results of the 2 solutions are consistent.

First of all, define all the provided formulas into rules and functions:

Clear["*"]
rule = ReleaseHold[
Hold[{f[r] =
Log[r - rp]/(2 KP) - Log[r - rN]/(2 KN) -
A Log[r^2 + r P + q] + ((2 (B + A P)) (ArcTan[(2 r + P)/Sqrt[4 q - P^2]] - π/
2))/Sqrt[4 q - P^2] - 20,
A = (R^2 (rp + rN) (rp^2 + rN^2 + 2 rp rN + R^2))/(
2 (3 rp^2 + rN^2 + 2 rp rN + R^2) (3 rN^2 + rp^2 + 2 rp rN + R^2)),
B = (R^2 (rp^2 + rN^2 + R^2) (rp^2 + rN^2 + rp rN + R^2))/((3 rp^2 + rN^2 +
2 rp rN + R^2) (3 rN^2 + rp^2 + 2 rp rN + R^2)),
KP = ((rp - rN) (3 rp^2 + rN^2 + 2 rp rN + R^2))/(2 R^2 rp^2),
KN = ((rp - rN) (3 rN^2 + rp^2 + 2 rp rN + R^2))/(2 R^2 rN^2),
rN = -(rp/3) - (2^(1/3) (3 rp^2 + 2 rp^4))/(
3 rp ((27 rp^2)/100 + 9 rp^4 + 7 rp^6 + Sqrt[
4 (3 rp^2 + 2 rp^4)^3 + ((27 rp^2)/100 + 9 rp^4 + 7 rp^6)^2])^(
1/3)) + ((27 rp^2)/100 + 9 rp^4 + 7 rp^6 + Sqrt[
4 (3 rp^2 + 2 rp^4)^3 + ((27 rp^2)/100 + 9 rp^4 + 7 rp^6)^2])^(1/3)/(
3 2^(1/3) rp), q = rp^2 + rN^2 + rp rN + R^2, P = rp + rN, rp = 4/10, R = 1}] /.
Set -> Rule];

f[r_] = Block[{f}, f[r] //. rule];

V[r_] = ((-(1/(50 r^3)) + 2 r + 1/(100 r^2 rp) + rp/r^2 + rp^3/r^2) (1 + 1/(100 r^2) +
r^2 - 1/(100 r rp) - rp/r - rp^3/r))/r /. rule;


Before moving on, it's worth analyzing the property of $f(r)$ and $v(r)$ first:

f[r] // N // Simplify
Limit[f[r], r -> 4/10]
Limit[f[r], r->Infinity]
Plot[f[x], {x, 4/10, 10}, PlotRange -> All]
(* -Infinity *)
(* -20 *)
(* -21.2263 + 0.780691 ArcTan[0.198677 + 0.942968 r] + 0.282187 Log[-0.4 + r] -
0.001025 Log[-0.0213856 + r] - 0.140581 Log[1.16901 + 0.421386 r + r^2] *)


It's not hard to notice $f(r)$ is a function defined in $(\frac{4}{10}, \infty)$ , monotone increasing from $-\infty$ to $-20$.

V[r] // N // Simplify
Plot[V[r], {r, -30, 20}]
(* (2. (0.01 - 0.489 r + 1. r^2 + r^4) (-0.01 + 0.2445 r + 1. r^4))/r^6 *)


$V(r)$ changes quickly near $0$.

Let's go on solving the PDE. Your equation is a first order PDE with a extreme variable coefficient, to solve it, one of the troublesome parts is determining $r$.

Since $f(r)$ is monotone and $f(r)=\frac{v-u}{2}$, we have $r=f^{(-1)}(\frac{v-u}{2})$. To calculate the inverse function of f, a direct use of InverseFunction is inefficient, so I choose interpolation for this task:

interinversef =
Interpolation[Transpose@{f@#, #} &@N[Range[4/10 + 10^-10, 100 + 10^-10, 1/100], 256]]


I've chosen a (probably unnecessary) high precision here to make sure the later calculation won't suffer from precision issue.

This is not the end, notice that I've only inversed the function from 4/10 + 10^-10 to 100 + 10^-10, because the variation of the function value near 4/10 is extremely fast and that after 100 is extremely slow and extents to infinity so interpolation is again insufficient. Here I choose 2 analytic approximation

frestpart1[r_] = f[r] /. Log[-(2/5) + r] -> 0;
fappro1[r_] = f[r] - frestpart1[r] + frestpart1[4/10];
inversefappro1[uv_] = InverseFunction[fappro1][uv] // First;

fappro2[r_] = f[r] /. Log[_] -> 0;
inversefappro2[uv_] = InverseFunction[fappro2][uv][[1]];


The selection of fappro2 isn't quite accurate actually, but I think it's enough and can't think out a better one currently.

Then we combine all these into one function:

With[{lowerbound = f[4/10 + 10^-10], upperbound = f[100 + 10^-10]},
inversef[uv_] :=
Piecewise[{{inversefappro2[uv], -20 > uv > upperbound}, {inversefappro1@uv,
uv < lowerbound}, {interinversef@uv, lowerbound <= uv <= upperbound}}]];


And solve the equation with NDSolve:

{u0 = 370, v0 = -10, uend = 400, vend = 20};
AbsoluteTiming[
sol = NDSolveValue[{(1/4) V[r] ψ[u, v] + Derivative[1, 1][ψ][u, v] ==
0, ψ[u, v0] == 0, ψ[u0, v] == Exp[-((v - vc)^2/(2 σ^2))]} /.
{vc -> 10, σ -> 3, r -> inversef[(v - u)/2]}, ψ, {u, u0, uend}, {v,
v0, vend}]]
Plot3D[sol[u, v], {u, u0, uend}, {v, v0, vend}, PlotRange -> All]


Hmm… no oscillation. Let's double-check it with FDM:

Clear@p
ugrid = 100; vgrid = 100; du = N[(uend - u0)/ugrid, 256]; dv = N[(vend - v0)/vgrid, 256];

formula = Subscript[ψ, E] - Subscript[ψ, S] + Subscript[ψ, W] -
1/32 V δu δv (-Subscript[u, E] - Subscript[u, N] + Subscript[v,
N] + Subscript[v, W]) (Subscript[ψ, E] + Subscript[ψ, W]);

vboundary[v_] = Exp[-((v - vc)^2/(2 σ^2))] /. {vc -> 10, σ -> N[3, 256]};

(p[{#, 0}] = N[0, 256]) & /@ Range[0, ugrid];
(p[{0, #}] = vboundary[v0 + # dv]) & /@ Range[0, vgrid];

With[{expr =
Function[{m, n}, #] &[
formula /.
Subscript[a_, b_] :> a[b] /.
{N -> {m + 1, n + 1}, S -> {m, n}, E -> {m, n + 1}, W -> {m + 1, n}} /.
{ψ -> p, δu -> du, δv -> dv, u -> First, v -> Last} /.
V -> V[inversef[(v0 + n dv - (u0 + m du))/2]] /.
{m -> m - 1, n -> n - 1}]},
p[{m_, n_}] := p[{m, n}] = expr[m, n]]

Block[{$MaxExtraPrecision = 500,$RecursionLimit = ∞}, p[{ugrid, vgrid}]];// AbsoluteTiming
dat = Table[p[{m, n}], {m, 0, ugrid}, {n, 0, vgrid}];

ListPlot3D[dat, PlotRange -> All]


Looks the same as the result from NDSolve.

# Epilog

I tend to guess my u0 and v0 isn't proper, currently value of (v-u)/2 is far from -20, while inversef varies quickly near -20.

To verify this guess, a smaller grid size is needed. The δu=δv=0.1 suggested by you is probably not enough given

inversef[-20.3] // N
inversef[-20.03] // N
inversef[-20.003] // N

(* 3.24399 *)
(* 33.3234 *)
(* 275.757 *)
`

I think the verification won't be too hard and I'll leave it to you. Now I'd like to stop here and go to bed.

• @xzczd...Thank you very much for your help.. – mat Nov 29 '15 at 9:32
• @mat You are welcome. If you feel satisfied with the answer, you can click the click the checkmark sign in the left-top corner of the answer to accept it. – xzczd Nov 29 '15 at 9:42