Revisions to Nonlinear system of ODEs with boundary conditions

logic error corrected

Source Link

edited Aug 11, 2018 at 18:15

62.1k
18
92
160

Error in boundary condition corrected

To solve this ODE system, begin by determining its behavior near r == 0.

In other words, g[r] vanishes to all orders at the origin for most values of n. However, if n == 1 is a positive integer, the same set of commands leaves g'[0]Derivative[n][g][0] undetermined. So, ng[eps] must be equalis approximately proportional to unityeps^n. Next, solve the equations numerically by shooting, using the method previously employed here and elsewhere on this site. (As noted in my comment above, DSolve returns unevaluated when applied to these equations.) A typical calculation for n == 5 is given by

eps = 10^-5; end = 20;
sp = ParametricNDSolveValue[{g'[r] == a[r] g[r]/r, a'[r]/r == g[r]^2 - 1, 
    a[eps] == 1n0, g[eps] == epseps^n0 gp, 
    WhenEvent[g[r] > 101/100, {bool = 1, "StopIntegration"}], 
    WhenEvent[a[r] < 0, {bool = 0, "StopIntegration"}]}, 
    {a[r], g[r]}, {r, eps, end + 5}, {gp, wp0, n0}, 
    WorkingPrecision -> wp0, Method -> "StiffnessSwitching"];
bl = 1/2;0; bu = 2; imax = 100; wp = 75; n = 5;
Do[bool = -1; bmiddle = (bl + bu)/2; st = sp[bmiddle, wp];wp, n]; 
    rm = (First[st] /. r -> "Domain")[[1, 2]]; 
    If[bool == 0, bl = bmiddle, bu = bmiddle]; ip = i; 
    If[bool == -1, Return[]], {i, imax}]; 
s[n] = st;

For larger values of n, wp must be increased. For instance, with n = 10, wp = 90 is required. Results for n = {1, 2, 5, 10} are plotted below.

Plot[{sts[1], s[2], s[5], s[10]}, {r, eps, Min[rm, end]10}, PlotRange -> All, AxesLabel -> {r, "a, g"}, 
    ImageSize -> Large, LabelStyle -> {Black, Bold, Medium}]

The solutions takes only a fewno more than several seconds each.

To solve this ODE system, begin by determining its behavior near r == 0.

In other words, g[r] vanishes to all orders at the origin for most values of n. However, if n == 1, the same set of commands leaves g'[0] undetermined. So, n must be equal to unity. Next, solve the equations numerically by shooting, using the method previously employed here and elsewhere on this site. (As noted in my comment above, DSolve returns unevaluated when applied to these equations.)

eps = 10^-5; end = 20;
sp = ParametricNDSolveValue[{g'[r] == a[r] g[r]/r, a'[r]/r == g[r]^2 - 1, 
    a[eps] == 1, g[eps] == eps gp, 
    WhenEvent[g[r] > 101/100, {bool = 1, "StopIntegration"}], 
    WhenEvent[a[r] < 0, {bool = 0, "StopIntegration"}]}, 
    {a[r], g[r]}, {r, eps, end + 5}, {gp, wp0}, 
    WorkingPrecision -> wp0, Method -> "StiffnessSwitching"];
bl = 1/2; bu = 2; imax = 100; wp = 75; 
Do[bool = -1; bmiddle = (bl + bu)/2; st = sp[bmiddle, wp]; 
    rm = (First[st] /. r -> "Domain")[[1, 2]]; 
    If[bool == 0, bl = bmiddle, bu = bmiddle]; ip = i; 
    If[bool == -1, Return[]], {i, imax}]; 
Plot[{st}, {r, eps, Min[rm, end]}, PlotRange -> All, AxesLabel -> {r, "a, g"}, 
    ImageSize -> Large, LabelStyle -> {Black, Bold, Medium}]

The solutions takes only a few seconds.

Error in boundary condition corrected

To solve this ODE system, begin by determining its behavior near r == 0.

In other words, g[r] vanishes to all orders at the origin for most values of n. However, if n is a positive integer, the same set of commands leaves Derivative[n][g][0] undetermined. So, g[eps] is approximately proportional to eps^n. Next, solve the equations numerically by shooting, using the method previously employed here and elsewhere on this site. (As noted in my comment above, DSolve returns unevaluated when applied to these equations.) A typical calculation for n == 5 is given by

eps = 10^-5; end = 20;
sp = ParametricNDSolveValue[{g'[r] == a[r] g[r]/r, a'[r]/r == g[r]^2 - 1, 
    a[eps] == n0, g[eps] == eps^n0 gp, 
    WhenEvent[g[r] > 101/100, {bool = 1, "StopIntegration"}], 
    WhenEvent[a[r] < 0, {bool = 0, "StopIntegration"}]}, 
    {a[r], g[r]}, {r, eps, end + 5}, {gp, wp0, n0}, 
    WorkingPrecision -> wp0, Method -> "StiffnessSwitching"];
bl = 0; bu = 2; imax = 100; wp = 75; n = 5;
Do[bool = -1; bmiddle = (bl + bu)/2; st = sp[bmiddle, wp, n]; 
    rm = (First[st] /. r -> "Domain")[[1, 2]]; 
    If[bool == 0, bl = bmiddle, bu = bmiddle]; ip = i; 
    If[bool == -1, Return[]], {i, imax}]; 
s[n] = st;

For larger values of n, wp must be increased. For instance, with n = 10, wp = 90 is required. Results for n = {1, 2, 5, 10} are plotted below.

Plot[{s[1], s[2], s[5], s[10]}, {r, eps, 10}, PlotRange -> All, AxesLabel -> {r, "a, g"}, 
    ImageSize -> Large, LabelStyle -> {Black, Bold, Medium}]

The solutions takes no more than several seconds each.

fixed minor bug in code

Source Link

edited Aug 9, 2018 at 13:02

bbgodfrey

62.1k
18
92
160

To solve this ODE system, begin by determining its behavior near r == 0.

Series[{g'[r] - (a[r] g[r])/r, a'[r]/r + 1 - g[r]^2}, {r, 0, 2}] // Normal;
Thread[Flatten@CoefficientList[r % /. a[0] -> n, r] == 0];
Solve[%, {a'[0], a''[0], a'''[0], a''''[0], g[0], g'[0], g''[0], g'''[0]}]

(* a'[0] -> 0, a''[0] -> -1, a'''[0] -> 0, a''''[0] -> 0,
   g[0] -> 0, g'[0] -> 0, g''[0] -> 0, g'''[0] -> 0  *)

In other words, g[r] vanishes to all orders at the origin for most values of n. However, if n == 1, the same set of commands leaves g'[0] undetermined. So, n must be equal to unity. Next, solve the equations numerically by shooting, using the method previously employed here and elsewhere on this site. (As noted in my comment above, DSolve returns unevaluated when applied to these equations.)

eps = 10^-5; end = 20;
sp = ParametricNDSolveValue[{g'[r] == a[r] g[r]/r, a'[r]/r == g[r]^2 - 1, 
    a[eps] == 1, g[eps] == eps gp, 
    WhenEvent[g[r] > 12101/10100, {bool = 1, "StopIntegration"}], 
    WhenEvent[{g[r] < 8/10, a[r]WhenEvent[a[r] < 0}, {bool = 0, "StopIntegration"}]}, 
    {a[r], g[r]}, {r, eps, end + 5}, {gp, wp0}, 
    WorkingPrecision -> wp0, Method -> "StiffnessSwitching"];
bl = 1/2; bu = 2; imax = 100; wp = 60;75; 
Do[bool = -1; bmiddle = (bl + bu)/2; st = sp[bmiddle, wp]; 
    rm = (First[st] /. r -> "Domain")[[1, 2]]; 
    If[bool == 0, bl = bmiddle, bu = bmiddle]; ip = i; 
    If[bool == -1, Return[]], {i, imax}]; 
Plot[{st}, {r, eps, Min[rm, end]}, PlotRange -> All, AxesLabel -> {r, "a, g"}, 
    ImageSize -> Large, LabelStyle -> {Black, Bold, Medium}]

The solutions takes only a few seconds.

To solve this ODE system, begin by determining its behavior near r == 0.

Series[{g'[r] - (a[r] g[r])/r, a'[r]/r + 1 - g[r]^2}, {r, 0, 2}] // Normal;
Thread[Flatten@CoefficientList[r % /. a[0] -> n, r] == 0];
Solve[%, {a'[0], a''[0], a'''[0], a''''[0], g[0], g'[0], g''[0], g'''[0]}]

(* a'[0] -> 0, a''[0] -> -1, a'''[0] -> 0, a''''[0] -> 0,
   g[0] -> 0, g'[0] -> 0, g''[0] -> 0, g'''[0] -> 0  *)

In other words, g[r] vanishes to all orders at the origin for most values of n. However, if n == 1, the same set of commands leaves g'[0] undetermined. So, n must be equal to unity. Next, solve the equations numerically by shooting, using the method previously employed here and elsewhere on this site. (As noted in my comment above, DSolve returns unevaluated when applied to these equations.)

eps = 10^-5; end = 20;
sp = ParametricNDSolveValue[{g'[r] == a[r] g[r]/r, a'[r]/r == g[r]^2 - 1, 
    a[eps] == 1, g[eps] == eps gp, 
    WhenEvent[g[r] > 12/10, {bool = 1, "StopIntegration"}], 
    WhenEvent[{g[r] < 8/10, a[r] < 0}, {bool = 0, "StopIntegration"}]}, 
    {a[r], g[r]}, {r, eps, end + 5}, {gp, wp0}, Method -> "StiffnessSwitching"];
bl = 1/2; bu = 2; imax = 100; wp = 60; 
Do[bool = -1; bmiddle = (bl + bu)/2; st = sp[bmiddle, wp]; 
    rm = (First[st] /. r -> "Domain")[[1, 2]]; 
    If[bool == 0, bl = bmiddle, bu = bmiddle]; ip = i; 
    If[bool == -1, Return[]], {i, imax}]; 
Plot[{st}, {r, eps, Min[rm, end]}, PlotRange -> All, AxesLabel -> {r, "a, g"}, 
    ImageSize -> Large, LabelStyle -> {Black, Bold, Medium}]

The solutions takes only a few seconds.

To solve this ODE system, begin by determining its behavior near r == 0.

Series[{g'[r] - (a[r] g[r])/r, a'[r]/r + 1 - g[r]^2}, {r, 0, 2}] // Normal;
Thread[Flatten@CoefficientList[r % /. a[0] -> n, r] == 0];
Solve[%, {a'[0], a''[0], a'''[0], a''''[0], g[0], g'[0], g''[0], g'''[0]}]

(* a'[0] -> 0, a''[0] -> -1, a'''[0] -> 0, a''''[0] -> 0,
   g[0] -> 0, g'[0] -> 0, g''[0] -> 0, g'''[0] -> 0  *)

In other words, g[r] vanishes to all orders at the origin for most values of n. However, if n == 1, the same set of commands leaves g'[0] undetermined. So, n must be equal to unity. Next, solve the equations numerically by shooting, using the method previously employed here and elsewhere on this site. (As noted in my comment above, DSolve returns unevaluated when applied to these equations.)

eps = 10^-5; end = 20;
sp = ParametricNDSolveValue[{g'[r] == a[r] g[r]/r, a'[r]/r == g[r]^2 - 1, 
    a[eps] == 1, g[eps] == eps gp, 
    WhenEvent[g[r] > 101/100, {bool = 1, "StopIntegration"}], 
    WhenEvent[a[r] < 0, {bool = 0, "StopIntegration"}]}, 
    {a[r], g[r]}, {r, eps, end + 5}, {gp, wp0}, 
    WorkingPrecision -> wp0, Method -> "StiffnessSwitching"];
bl = 1/2; bu = 2; imax = 100; wp = 75; 
Do[bool = -1; bmiddle = (bl + bu)/2; st = sp[bmiddle, wp]; 
    rm = (First[st] /. r -> "Domain")[[1, 2]]; 
    If[bool == 0, bl = bmiddle, bu = bmiddle]; ip = i; 
    If[bool == -1, Return[]], {i, imax}]; 
Plot[{st}, {r, eps, Min[rm, end]}, PlotRange -> All, AxesLabel -> {r, "a, g"}, 
    ImageSize -> Large, LabelStyle -> {Black, Bold, Medium}]

The solutions takes only a few seconds.

added missing result

Source Link

edited Aug 8, 2018 at 23:12

bbgodfrey

62.1k
18
92
160

To solve this ODE system, begin by determining its behavior near r == 0.

Series[{g'[r] - (a[r] g[r])/r, a'[r]/r + 1 - g[r]^2}, {r, 0, 2}] // Normal;
Thread[Flatten@CoefficientList[r % /. a[0] -> n, r] == 0];
Solve[%, {a'[0], a''[0], a'''[0], a''''[0], g[0], g'[0], g''[0], g'''[0]}] 

(* a'[0] -> 0, a''[0] -> -1, a'''[0] -> 0, a''''[0] -> 0,
   g[0] -> 0, g'[0] -> 0, g''[0] -> 0, g'''[0] -> 0  *)

In other words, g[r] must vanishvanishes to all orders at the origin for most values of n. However, if n == 1, the same set of commands leaves g'[0] undetermined. So, n must be equal to unity. Next, solve the equations numerically by shooting, using the method previously employed here and elsewhere on this site. (As noted in my comment above, DSolve returns unevaluated when applied to these equations.)

eps = 10^-5; end = 20;
sp = ParametricNDSolveValue[{g'[r] == a[r] g[r]/r, a'[r]/r == g[r]^2 - 1, 
    a[eps] == 1, g[eps] == eps gp, 
    WhenEvent[g[r] > 12/10, {bool = 1, "StopIntegration"}], 
    WhenEvent[{g[r] < 8/10, a[r] < 0}, {bool = 0, "StopIntegration"}]}, 
    {a[r], g[r]}, {r, eps, end + 5}, {gp, wp0}, Method -> "StiffnessSwitching"];
bl = 1/2; bu = 2; imax = 100; wp = 60; 
Do[bool = -1; bmiddle = (bl + bu)/2; st = sp[bmiddle, wp]; 
    rm = (First[st] /. r -> "Domain")[[1, 2]]; 
    If[bool == 0, bl = bmiddle, bu = bmiddle]; ip = i; 
    If[bool == -1, Return[]], {i, imax}]; 
Plot[{st}, {r, eps, Min[rm, end]}, PlotRange -> All, AxesLabel -> {r, "a, g"}, 
    ImageSize -> Large, LabelStyle -> {Black, Bold, Medium}]

The solutions takes only a few seconds.

To solve this ODE system, begin by determining its behavior near r == 0.

Series[{g'[r] - (a[r] g[r])/r, a'[r]/r + 1 - g[r]^2}, {r, 0, 2}] // Normal;
Thread[Flatten@CoefficientList[r % /. a[0] -> n, r] == 0];
Solve[%, {a'[0], a''[0], a'''[0], a''''[0], g[0], g'[0], g''[0], g'''[0]}]

In other words, g[r] must vanish to all orders at the origin for most values of n. However, if n == 1, the same set of commands leaves g'[0] undetermined. So, n must be equal to unity. Next, solve the equations numerically by shooting, using the method previously employed here and elsewhere on this site.

eps = 10^-5; end = 20;
sp = ParametricNDSolveValue[{g'[r] == a[r] g[r]/r, a'[r]/r == g[r]^2 - 1, 
    a[eps] == 1, g[eps] == eps gp, 
    WhenEvent[g[r] > 12/10, {bool = 1, "StopIntegration"}], 
    WhenEvent[{g[r] < 8/10, a[r] < 0}, {bool = 0, "StopIntegration"}]}, 
    {a[r], g[r]}, {r, eps, end + 5}, {gp, wp0}, Method -> "StiffnessSwitching"];
bl = 1/2; bu = 2; imax = 100; wp = 60; 
Do[bool = -1; bmiddle = (bl + bu)/2; st = sp[bmiddle, wp]; 
    rm = (First[st] /. r -> "Domain")[[1, 2]]; 
    If[bool == 0, bl = bmiddle, bu = bmiddle]; ip = i; 
    If[bool == -1, Return[]], {i, imax}]; 
Plot[{st}, {r, eps, Min[rm, end]}, PlotRange -> All, AxesLabel -> {r, "a, g"}, 
    ImageSize -> Large, LabelStyle -> {Black, Bold, Medium}]

The solutions takes only a few seconds.

To solve this ODE system, begin by determining its behavior near r == 0.

Series[{g'[r] - (a[r] g[r])/r, a'[r]/r + 1 - g[r]^2}, {r, 0, 2}] // Normal;
Thread[Flatten@CoefficientList[r % /. a[0] -> n, r] == 0];
Solve[%, {a'[0], a''[0], a'''[0], a''''[0], g[0], g'[0], g''[0], g'''[0]}] 

(* a'[0] -> 0, a''[0] -> -1, a'''[0] -> 0, a''''[0] -> 0,
   g[0] -> 0, g'[0] -> 0, g''[0] -> 0, g'''[0] -> 0  *)

In other words, g[r] vanishes to all orders at the origin for most values of n. However, if n == 1, the same set of commands leaves g'[0] undetermined. So, n must be equal to unity. Next, solve the equations numerically by shooting, using the method previously employed here and elsewhere on this site. (As noted in my comment above, DSolve returns unevaluated when applied to these equations.)

eps = 10^-5; end = 20;
sp = ParametricNDSolveValue[{g'[r] == a[r] g[r]/r, a'[r]/r == g[r]^2 - 1, 
    a[eps] == 1, g[eps] == eps gp, 
    WhenEvent[g[r] > 12/10, {bool = 1, "StopIntegration"}], 
    WhenEvent[{g[r] < 8/10, a[r] < 0}, {bool = 0, "StopIntegration"}]}, 
    {a[r], g[r]}, {r, eps, end + 5}, {gp, wp0}, Method -> "StiffnessSwitching"];
bl = 1/2; bu = 2; imax = 100; wp = 60; 
Do[bool = -1; bmiddle = (bl + bu)/2; st = sp[bmiddle, wp]; 
    rm = (First[st] /. r -> "Domain")[[1, 2]]; 
    If[bool == 0, bl = bmiddle, bu = bmiddle]; ip = i; 
    If[bool == -1, Return[]], {i, imax}]; 
Plot[{st}, {r, eps, Min[rm, end]}, PlotRange -> All, AxesLabel -> {r, "a, g"}, 
    ImageSize -> Large, LabelStyle -> {Black, Bold, Medium}]

The solutions takes only a few seconds.

Source Link

created Aug 8, 2018 at 23:05

bbgodfrey

62.1k
18
92
160

Loading

Stack Exchange Network

Return to Answer