Is there a way to automate this procedure?

Question

I have some code which numerically solves a differential equation for a given value of the parameter $M0$ and for some initial conditions. The initial conditions will remain fixed throughout the procedure and we want to vary the value of the parameter $M0$.

The code is given below

eqnofmotion := 1/r^3 D[r^3 D[f[r], r], r] + M0^2/(r^2 + L^2)^2 f[r]
L := 1
M0 := 2 Sqrt[2]
sltn = NDSolve[{eqnofmotion == 0, f[10^-5] == 1, f'[10^-5] == 0}, 
   f, {r, 10^-5, 10^2}];
Plot[ f[r] /. sltn, {r, 10^-5, 10^2}, PlotRange -> {-1, 1}, 
 PlotStyle -> {Thick, Black}, 
 BaseStyle -> {18, FontFamily -> "Times New Roman"}, 
 AxesLabel -> {"\[Rho]", "\[Delta]L(\[Rho])"}, 
 PlotLegends -> Automatic]
f[r] /. sltn /. r :> 10^2

After this, the following procedure takes place. Define as a function the solution above, and its derivative and evaluate them at $100$.

function0[r_] := f[r] /. sltn
derivative0[r_] := function0'[r]
function0[r] /. r :> 100
derivative0[r] /. r :> 100

From the above values, we can calculate the operator and the source of the state. The code is given below

scalarasymptote[r_] := operator/r^3 + source r
(*Define and solve the system to determine operator and source*)
\
Solve[{scalarasymptote[100] == 0.0000999607, 
  scalarasymptote'[100] == - 1.99889 10^-6}, {operator, source}]

And now that we have the values of the operator and the source, we can calculate the square of the coupling, which is the $a$ in the following code,

(*Now use that \[ScriptCapitalJ] = g^2/\[CapitalLambda]^2 \
\[ScriptCapitalO]*) 
(*\[CapitalLambda] is the UV cutoff.*)
\[CapitalLambda] := 10^2
NSolve[ 2.49983 10^-7 == 74.9624 a/\[CapitalLambda]^2, a]

The goal is to repeat all of the above steps for many values of $M0$, that are smaller than $2 \sqrt{2}$. Create a list of the square of all the values of $M0$ and the corresponding $a$ values and plot the data.

The way I do it, is to copy and paste the above piece of code, and run it each time for a different value of $M0$, for which I use another symbol to save it, find $a$ each time, until I reach the value $M0=0$ and then plot the data.

The thing I would like to know is if and how I can automate the above procedure or at least some steps of it, in order to get more points a bit more easily and within a reasonable amount of time. Copying and pasting and running the whole thing each time takes quite some time and the goal is to do this for excited states as well, start with higher $M0$ and go all the way down to $0$.

Thanks in advance.

Answer 1

Try writing your code as a series of functions that link back to each other rather than manually carrying results through. For example, instead of copying 0.0000999607 from f[100]'s result, just write f[100].

Some stuff does not need to be redefined every time it appears. I find it's usually best to separate these out:

eqnofmotion[L_, M0_] := 
  1/r^3 D[r^3 D[f[r], r], r] + M0^2/(r^2 + L^2)^2 f[r];
scalarasymptote[r_] := operator/r^3 + source r;

I have rewritten eqnofmotion to explicitly depend on L and M0, so as to facilitate substituting these values in later.

sltn[L_, M0_] := 
  NDSolve[{eqnofmotion[L, M0] == 0, f[10^-5] == 1, f'[10^-5] == 0}, 
   f, {r, 10^-5, 10^2}];

The first sltn can be written so that you can pass L and M0 into it. This requires both it and eqnofmotion to depend on L and M0.

sltn2[L_, M0_] := 
  Solve[Evaluate[{scalarasymptote[100] == f[100], 
      scalarasymptote'[100] == f'[100]} /. 
     First@sltn[L, M0]], {operator, source}];

Then we write sltn2, also dependent on L and M0, so that it automatically incorporates the values found from the first solution with /. First@sltn[L, M0].

sltn3[L_, M0_, Λ_] := 
  Solve[(source == operator a/Λ^2) /. 
    First@sltn2[L, M0], a];

sltn3 works much the same way, though it would probably be better to manually evaluate the solution first and then substitute the values in in this case. That would avoid calling Solve unnecessarily, and the general solver functions are typically a bit slower. However, in this case, it doesn't really matter.

With this, we have a function that directly relates L, M0, and Λ to a. In fact, we can plot it directly in Mathematica:

Plot[a /. sltn3[1, m, 100], {m, 0, 2 Sqrt[2]}]

a starts fairly level when M0 equals 0 and drops off sharply near 2 Sqrt[2].