Partial differential equation, Finite difference Method

Question

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.

Please note that here

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.

Answer 1

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.