0
$\begingroup$

I have a code given by

V =((1/2 - (3 r)/2 + r^2) (r^8 (-1 + (3 r)/2) + 
4 (2 - (9 r)/2 + 2 r^2) + 2 r^4 (3 - 9 r + 6 r^2)))/(r^6 (2 + 
r^4)^2);

X= r + 5/2 ArcTan[4 (-(3/4) + r)] + 3/4 Log[1/2 - (3 r)/2 + r^2];

zz = X /. r -> 3/2

fV[z0_?NumericQ] :=With[{z = SetPrecision[z0 + zz, 100 + 1]},If[Abs[z] 
<= 35,Re[V /. FindRoot[X == z, {r, 10000001/10000000}, 
  MaxIterations -> 10000, WorkingPrecision -> 100]], 0]]

My problem: consider this formula*

Subscript[\[Psi], m, n, p] = 
Subscript[\[Psi], m - 1, n - 1, p] + Subscript[\[Psi], m, n, p - 1] -
Subscript[\[Psi], m - 1, n - 1, p - 1] - 
h^2/8 (fV[(p - (n - 1))/2 h] Subscript[\[Psi], m - 1, n - 1, p] + 
 fV[((p - 1) - n)/2 h] Subscript[\[Psi], m, n, p - 1])

with the following initial conditions

Subscript[\[Psi], 0, 0, 0] = 1;
Subscript[\[Psi], 1, 1, 0] = 1;
Subscript[\[Psi], 2, 2, 0] = 1;
Subscript[\[Psi], 3, 3, 0] = 1;
Subscript[\[Psi], 4, 4, 0] = 1;
Subscript[\[Psi], 5, 5, 0] = 1;
Subscript[\[Psi], 6, 6, 0] = 1;
Subscript[\[Psi], 7, 7, 0] = 1;
Subscript[\[Psi], 8, 8, 0] = 1;
Subscript[\[Psi], 9, 9, 0] = 1;
Subscript[\[Psi], 10, 10, 0] = 1;
Subscript[\[Psi], 0, 0, p] = Exp[-0.25 (p h)^2];
h = .5;

I need a table given by

Table[{t, Subscript[\[Psi], t/h, t/h, t/h]}, {t, 0, 5, h}]

as the result of calculations. How can I tell the Mathematica to calculate Subscript[\[Psi], m, n, p] from formula* and provide me the mentioned table? I will be thankful if someone helps.

$\endgroup$
3
$\begingroup$

You could define a function instead of using Subscript (this is in general a good idea) and employ memoization (in order to preventing that values have to be computed multiple times:

ClearAll[fV];
fV[z0_?NumericQ] := fV[z0] = With[{z = SetPrecision[z0 + zz, 100 + 1]}, 
   If[Abs[z] <= 35, 
    Re[V /. FindRoot[X == z, {r, 10000001/10000000}, 
       MaxIterations -> 10000, WorkingPrecision -> 100]], 0]];

h = 1/2;
ClearAll[Ψ];
Ψ[0, 0, 0] = 1;
Ψ[1, 1, 0] = 1;
Ψ[2, 2, 0] = 1;
Ψ[3, 3, 0] = 1;
Ψ[4, 4, 0] = 1;
Ψ[5, 5, 0] = 1;
Ψ[6, 6, 0] = 1;
Ψ[7, 7, 0] = 1;
Ψ[8, 8, 0] = 1;
Ψ[9, 9, 0] = 1;
Ψ[10, 10, 0] = 1;
Ψ[0, 0, p_] := Ψ[0, 0, p] = Exp[- 1/4 (p h)^2];
Ψ[m_, n_, p_] := Ψ[m, n, p] = Ψ[m - 1, n - 1, p] + Ψ[m, n,  - 1] - Ψ[m - 1, n - 1, p - 1] - h^2/8 (fV[(p - (n - 1))/2 h] Ψ[m - 1, n - 1, p] + fV[((p - 1) - n)/2 h] Ψ[m, n, p - 1])
Table[{t, Ψ[t/h, t/h, t/h]}, {t, 0, 5, h}]

{{0, 1}, {1/2, 0.936178}, {1, 0.766748}, {3/2, 0.54557}, {2, 0.330802}, {5/2, 0.161054}, {3, 0.0478391}, {7/2, -0.0171044}, {4, -0.0494697}, {9/2, -0.0633942}, {5, -0.068163}}

$\endgroup$
  • $\begingroup$ Thank you very much. When I employ it for higher values of t, say t=21, with h=1/10 and \[CapitalPsi][m_, n_, 0] := \[CapitalPsi][m, n, 0] = 1;, I cannot make its table and it gives some FindRoot errors. Could you please tell me why? $\endgroup$ – Mehrab Feb 22 '19 at 13:56
  • $\begingroup$ Consider h=1/50 $\endgroup$ – Mehrab Feb 22 '19 at 14:30
  • $\begingroup$ I obtain plenty of recursion limit errors. Probably your initial conditions are not sufficient; the recursive evaluation of Ψ does not end up in one of the cases Ψ[0, 0, 0], Ψ[1, 1, 0],...., or Ψ[0, 0, p]. To Ψ calls itseld recursively over and over again infinitely often. Of course, a computer cannot do that so Mathematica breaks this vicious circle at some point and throws an error. $\endgroup$ – Henrik Schumacher Feb 22 '19 at 15:19

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.