Skip to main content
added 235 characters in body
Source Link
Mariusz Iwaniuk
  • 16.5k
  • 1
  • 27
  • 46

Another way:

eps = 10^-10; end = 10;
sol[n_] := NDSolve[
{g'[r] == a[r]*g[r]/r, a'[r]/r == g[r]^2 - 1, a[eps] == n, g[end] == 1}
,{a[r], g[r]}, {r, eps, end}, Method -> {"Shooting", 
"StartingInitialConditions" -> {a[end] == eps, g[end] == (1 - eps)}}];

Plot[Evaluate[g[r] /. sol[#] & /@ {1, 2, 5, 10}], {r, eps, end},PlotRange -> All, 
PlotLegends -> {"n=1", "n=2", "n=5", "n=10"}, PlotLabels -> "g(r)"]
Plot[Evaluate[a[r] /. sol[#] & /@ {1, 2, 5, 10}], {r, eps, end}, 
 PlotRange -> All,  
PlotLegends -> {"n=1", "n=2", "n=5", "n=10"}, 
 PlotLabels -> "a(r)"]

enter image description here

For larger values e.g.end=20:

eps = 10^-20; end = 20;
sol[n_] := NDSolve["Code the same",WorkingPrecision -> 25];
Plot["Code the same"] // Quiet(*Error messages form InterpolatingFunction to Quiet *)

Another way:

eps = 10^-10; end = 10;
sol[n_] := NDSolve[
{g'[r] == a[r]*g[r]/r, a'[r]/r == g[r]^2 - 1, a[eps] == n, g[end] == 1}
,{a[r], g[r]}, {r, eps, end}, Method -> {"Shooting", 
"StartingInitialConditions" -> {a[end] == eps, g[end] == (1 - eps)}}];

Plot[Evaluate[g[r] /. sol[#] & /@ {1, 2, 5, 10}], {r, eps, end}, 
PlotLegends -> {"n=1", "n=2", "n=5", "n=10"}, PlotLabels -> "g(r)"]
Plot[Evaluate[a[r] /. sol[#] & /@ {1, 2, 5, 10}], {r, eps, end}, 
 PlotRange -> All, PlotLegends -> {"n=1", "n=2", "n=5", "n=10"}, 
 PlotLabels -> "a(r)"]

enter image description here

Another way:

eps = 10^-10; end = 10;
sol[n_] := NDSolve[
{g'[r] == a[r]*g[r]/r, a'[r]/r == g[r]^2 - 1, a[eps] == n, g[end] == 1}
,{a[r], g[r]}, {r, eps, end}, Method -> {"Shooting", 
"StartingInitialConditions" -> {a[end] == eps, g[end] == (1 - eps)}}];

Plot[Evaluate[g[r] /. sol[#] & /@ {1, 2, 5, 10}], {r, eps, end},PlotRange -> All, 
PlotLegends -> {"n=1", "n=2", "n=5", "n=10"}, PlotLabels -> "g(r)"]
Plot[Evaluate[a[r] /. sol[#] & /@ {1, 2, 5, 10}], {r, eps, end}, PlotRange -> All,  
PlotLegends -> {"n=1", "n=2", "n=5", "n=10"}, PlotLabels -> "a(r)"]

enter image description here

For larger values e.g.end=20:

eps = 10^-20; end = 20;
sol[n_] := NDSolve["Code the same",WorkingPrecision -> 25];
Plot["Code the same"] // Quiet(*Error messages form InterpolatingFunction to Quiet *)
Source Link
Mariusz Iwaniuk
  • 16.5k
  • 1
  • 27
  • 46

Another way:

eps = 10^-10; end = 10;
sol[n_] := NDSolve[
{g'[r] == a[r]*g[r]/r, a'[r]/r == g[r]^2 - 1, a[eps] == n, g[end] == 1}
,{a[r], g[r]}, {r, eps, end}, Method -> {"Shooting", 
"StartingInitialConditions" -> {a[end] == eps, g[end] == (1 - eps)}}];

Plot[Evaluate[g[r] /. sol[#] & /@ {1, 2, 5, 10}], {r, eps, end}, 
PlotLegends -> {"n=1", "n=2", "n=5", "n=10"}, PlotLabels -> "g(r)"]
Plot[Evaluate[a[r] /. sol[#] & /@ {1, 2, 5, 10}], {r, eps, end}, 
PlotRange -> All, PlotLegends -> {"n=1", "n=2", "n=5", "n=10"}, 
PlotLabels -> "a(r)"]

enter image description here