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.
ψfrom it. $\endgroup$Rremains undefined, and the question now contains two inconsistent definitions forrN. Yu do not make it easy for people to help you. $\endgroup$u0andv0, the image I obtained is quite different from yours. $\endgroup$