Skip to main content

differential equations - Solution of nonlinear system with boundary conditions


I'm try solve the following coupled ODEs with boundary conditions:


$I. \ \ \ \ \ \ \ \dfrac{4}{r}[1+A(r)]\left(\dfrac{dH}{dr}\right)^2+\dfrac{dG}{dr}=0$


$II. \ \ \ \ \dfrac{1}{r}\left(\dfrac{dA}{dr}\right)+F(r)+k^2G(r)=0$



$III. \ \ \ \ \left(\dfrac{dG}{dr}\right)^2+4k^2\left(G(r)+\dfrac{F(r)}{k^2}\right)^2\left(\dfrac{dH}{dr}\right)^2-1.6H(r)\left(\dfrac{dH}{dr}\right)^2=0$


$IV. \ \ \ \ \dfrac{d^2F}{dr^2}+\dfrac{1}{r}\dfrac{dF}{dr}+\dfrac{1}{r}\dfrac{dA}{dr}-4k^2F(r)\left(\dfrac{dH}{dr}\right)^2=0$


where $k>0$. The boundary conditions are


$H(0)=1\\A(0)=0\\G(L)=0\\F'(0)=0, \ \ F(L)=0$


with $L$ being the boundary, such that can be adjusted in order to satisfy the boundary conditions above for each value of parameter $k$ choosed. For this physical problem, it is desirable that the ODEs of first order obey other conditions:


$H'(0)=0, \ \ H(L)=H'(L)=0\\A'(0)=0, \ \ A'(L)=0$


Thus, from ODEs, it is possible to conclude


The ODE $I$ guarantees $G'(0)=G'(L)=0$


The ODE $II$ guarantees $F(L)=0$


I tried solve, for example, with k=1 using the shooting method. In this case, it was only possible to execute up to $L = 1.45$ (end=1.45).



k = 1; eps = 10^-10; end = 1.45;
ode1 = 4 (h'[r] ^2) (a[r] + 1)/r == -G'[r];
ode2 = a'[r]/r == -(k^2) (G[r]) - f[r];
ode3 = (G'[r])^2 + 4 (k^2) (h'[r]^2) (G[r] + f[r]/k^2)^2 == (1.6) h[r] (h'[r]^2);
ode4 = f''[r] + f'[r]/r + (a'[r]/r) == 4 (k^2) (f[r]) ((h'[r])^2);
bcs = {h[eps] == (1 - eps), a[eps] == eps, G[end] == eps, f'[eps] == eps, f[end] == eps};

sol = NDSolve[{ode1, ode2, ode3, ode4, bcs}, {a, h, f, G}, {r, eps, end},
Method -> {"Shooting","StartingInitialConditions" -> {f[end] == eps, G[end] == eps}}];


Plot[Evaluate[{h[r], h'[r]} /. sol[[3]]], {r, eps, end},
PlotLegends -> {"H", "H'"}, PlotRange -> All]
Plot[Evaluate[{a[r], a'[r]} /. sol[[3]]], {r, eps, end},
PlotLegends -> {"A", "A'"}, PlotRange -> All]
Plot[Evaluate[{G[r], G'[r]} /. sol[[3]]], {r, eps, end},
PlotLegends -> {"G", "G'"}, PlotRange -> All]
Plot[Evaluate[{f[r], f'[r]} /. sol[[3]]], {r, eps, end},
PlotLegends -> {"F", "F'"}, PlotRange -> All]

enter image description here



However, for values $L>1.45$ appear the error message



"NDSolve::ndsz: At r == 1.294087660532624`, step size is effectively zero; singularity or stiff system suspected.."



But I believe it is possible to find the correct value of $L$ that satisfies all boundary conditions. The behavior of the functions by the plots suggests that this value can be near to 2 or 3. I can not resolve this error message for this problem.


Added_1: Specifically, I am interested in obtaining numerical solutions for the cases of $k=0.1$, $k=1$ and $k=2$. Then, there must be a value $L$ for each $k$ that satisfies all the boundary conditions of the problem.


Any help or altenative solution is welcome.



Answer



There is a solution satisfying all the boundary conditions for L = 2.665. I'll show you how to build this solution. First, we solve explicitly the first and third equations for the derivatives, we have


ode1 = 4 (h'[r]^2) (a[r] + 1)/r == -G'[r];

ode3 = (G'[r])^2 + 
    4 (k^2) (h'[r]^2) (G[r] + f[r]/k^2)^2 == (8/5) h[r] (h'[r]^2);
s = Solve[{ode1, ode3}, {h'[r], G'[r]}] // FullSimplify

Here we have three pairs of roots:


{{Derivative[1][h][r] -> 0, 
  Derivative[1][G][r] -> 
   0}, {Derivative[1][h][r] -> -((
    I r Sqrt[5 (f[r] + k^2 G[r])^2 - 2 k^2 h[r]])/(
    2 Sqrt[5] Sqrt[k^2 (1 + a[r])^2])), 

  Derivative[1][G][r] -> (5 r (f[r] + k^2 G[r])^2 - 2 k^2 r h[r])/(
   5 k^2 (1 + a[r]))}, {Derivative[1][h][r] -> (
   I r Sqrt[5 (f[r] + k^2 G[r])^2 - 2 k^2 h[r]])/(
   2 Sqrt[5] Sqrt[k^2 (1 + a[r])^2]), 
  Derivative[1][G][r] -> (5 r (f[r] + k^2 G[r])^2 - 2 k^2 r h[r])/(
   5 k^2 (1 + a[r]))}}

The first pair of roots corresponds to the trivial solution $H(r)=1,G(r)=0$. On this branch the second and fourth equations become linear, and their solution can be obtained in an explicit form. We are interested in the second pair of roots, from which we will form a new system of equations:


ode13 = {Derivative[1][G][r] == (
    5 r (f[r] + k^2 G[r])^2 - 2 k^2 r h[r])/(5 k^2 (1 + a[r])), 

   Derivative[1][h][r] == -(( 
     r Sqrt[-5 (f[r] - k^2 G[r])^2 + 2 k^2 h[r]])/(
     2 Sqrt[5] Sqrt[k^2 (1 + a[r])^2]))};
ode24 = {a'[r]/r == -(k^2) (G[r]) - f[r], 
  f''[r] + f'[r]/r + (a'[r]/r) == 4 (k^2) (f[r]) ((h'[r])^2)};

We reduce this system to the dimensionless form by making a substitution r->L*r. Then the parameter $L$ becomes a parameter of the equations, and the system and boundary conditions have the form:


k = 1; eps = 10^-10; L = 1.638;
ode13 = {Derivative[1][G][r] == 
    L^2*(5 r (f[r] + k^2 G[r])^2 - 2 k^2 r h[r])/(5 k^2 (1 + a[r])), 

   Derivative[1][h][r] == -L^2*( 
     r Sqrt[-5 (f[r] + k^2 G[r])^2 + 2 k^2 h[r]])/(
     2 Sqrt[5] Sqrt[k^2 (1 + a[r])^2])};
ode24 = {a'[r]/r == -L^2*((k^2) (G[r]) - f[r]), 
    f''[r] + f'[r]/r + (a'[r]/r) == 4 (k^2) (f[r]) ((h'[r])^2)};
bcs = {h[eps] == 1, a[eps] == 0, G[1] == 0, f'[eps] == 0, f[1] == 0};
sol = NDSolve[{ode13, ode24, bcs}, {a, h, f, G}, {r, eps, 1}, 
   Method -> {"Shooting", 
     "StartingInitialConditions" -> {f[eps] == -0.1, G[eps] == 0.33}}];
{Plot[Evaluate[Re[{h[r], h'[r]}] /. sol], {r, eps, 1}, 

  PlotLegends -> {"H", "H'"}, PlotRange -> All],
 Plot[Evaluate[Re[{a[r], a'[r]}] /. sol], {r, eps, 1}, 
  PlotLegends -> {"A", "A'"}, PlotRange -> All],
 Plot[Evaluate[Re[{G[r], G'[r]}] /. sol], {r, eps, 1}, 
  PlotLegends -> {"G", "G'"}, PlotRange -> All],
 Plot[Evaluate[Re[{f[r], f'[r]}] /. sol], {r, eps, 1}, 
  PlotLegends -> {"F", "F'"}, PlotRange -> All]}

fig1


Another identity transformation of the system allows us to find a solution with L = 2:



k = 1; eps = 10^-10; L = 2;
ode = {G'[r] == (L^2 r (5 (f[r] + k^2 G[r])^2 - 2 k^2 h[r]))/(
    5 k^2 (1 + a[r])), 
   Derivative[1][h][r] == -((
     L^2 r Sqrt[-5 f[r]^2 - 10 k^2 f[r] G[r] - 5 k^4 G[r]^2 + 
       2 k^2 h[r]])/Sqrt[20 k^2 + 40 k^2 a[r] + 20 k^2 a[r]^2]), 
   r*f''[r] + f'[r] - (r*L^2*((k^2) (G[r]) + f[r])) == 
    r*4 (k^2) (f[r]) ((h'[r])^2), 
   a'[r] == -r*L^2*((k^2) (G[r]) + f[r])};
bcs = {h[eps] == 1, a[eps] == 0, G[1] == 0, f'[eps] == 0, f[1] == 0};


sol = NDSolve[{ode, bcs}, {a, h, f, G}, {r, eps, 1}, 
   Method -> {"Shooting", 
     "StartingInitialConditions" -> {f[eps] == -0.1, G[eps] == 0.33}}];


{Plot[Evaluate[Re[{h[r], h'[r]}] /. sol], {r, eps, 1}, 
  PlotLegends -> {"H", "H'"}, PlotRange -> All],
 Plot[Evaluate[Re[{a[r], a'[r]}] /. sol], {r, eps, 1}, 
  PlotLegends -> {"A", "A'"}, PlotRange -> All],

 Plot[Evaluate[Re[{G[r], G'[r]}] /. sol], {r, eps, 1}, 
  PlotLegends -> {"G", "G'"}, PlotRange -> All],
 Plot[Evaluate[Re[{f[r], f'[r]}] /. sol], {r, eps, 1}, 
  PlotLegends -> {"F", "F'"}, PlotRange -> All]}

fig2


Finally, we point out a solution with h[L] tending to zero. In this case L = 2.665


k = 1; eps = 10^-10; L = 2.665;
ode = {G'[r] == (L^2 r (5 (f[r] + k^2 G[r])^2 - 2 k^2 h[r]))/(
    5 k^2 (1 + a[r])), 

   Derivative[1][h][r] == -((
     L^2 r Sqrt[-5 f[r]^2 - 10 k^2 f[r] G[r] - 5 k^4 G[r]^2 + 
       2 k^2 h[r]])/Sqrt[20 k^2 + 40 k^2 a[r] + 20 k^2 a[r]^2]), 
   r*f''[r] + f'[r] - (r*L^2*((k^2) (G[r]) + f[r])) == 
    r*4 (k^2) (f[r]) ((h'[r])^2), 
   a'[r] == -r*L^2*((k^2) (G[r]) + f[r])};
bcs = {h[eps] == 1, a[eps] == 0, G[1] == 0, f'[eps] == 0, f[1] == 0};

sol = NDSolve[{ode, bcs}, {a, h, f, G}, {r, eps, 1}, 
   Method -> {"Shooting", 

     "StartingInitialConditions" -> {f[eps] == -0.340071, 
       G[eps] == 0.736822}}];


{Plot[Evaluate[Re[{h[r], h'[r]}] /. sol], {r, eps, 1}, 
  PlotLegends -> {"H", "H'"}, PlotRange -> All],
 Plot[Evaluate[Re[{a[r], a'[r]}] /. sol], {r, eps, 1}, 
  PlotLegends -> {"A", "A'"}, PlotRange -> All],
 Plot[Evaluate[Re[{G[r], G'[r]}] /. sol], {r, eps, 1}, 
  PlotLegends -> {"G", "G'"}, PlotRange -> All],

 Plot[Evaluate[Re[{f[r], f'[r]}] /. sol], {r, eps, 1}, 
  PlotLegends -> {"F", "F'"}, PlotRange -> All]}

fig3


We indicate the solution algorithm in the case k = 2. We set A = 1, we define two functions


k = 2; eps = 10^-10; L = 1;
ode = {G'[
     r] == (L^2 r (5 (f[r] + k^2 G[r])^2 - 2 k^2 h[r]))/(5 k^2 (1 + 
         a[r])), Derivative[1][h][
     r] == -((L^2 r Sqrt[-5 f[r]^2 - 10 k^2 f[r] G[r] - 

           5 k^4 G[r]^2 + 2 k^2 h[r]])/
       Sqrt[20 k^2 + 40 k^2 a[r] + 20 k^2 a[r]^2]), 
   r*f''[r] + f'[r] - (r*L^2*((k^2) (G[r]) + f[r])) == 
    r*4 (k^2) (f[r]) ((h'[r])^2), 
   a'[r] == -r*L^2*((k^2) (G[r]) + f[r])};
 bcsp = {h[eps] == 1, a[eps] == 0, f'[eps] == 0, f[eps] == p, 
   G[eps] == q};
g = ParametricNDSolveValue[{ode, bcsp}, G, {r, eps, 1}, {p, q}, 
  Method -> {"StiffnessSwitching", "NonstiffTest" -> False}, 
  MaxSteps -> Infinity]


 F = 
 ParametricNDSolveValue[{ode, bcsp}, f, {r, eps, 1}, {p, q}, 
  Method -> {"StiffnessSwitching", "NonstiffTest" -> False}, 
  MaxSteps -> Infinity]

To fulfill the boundary conditions for r = 1, we find the roots of the system of equations


FindRoot[{g[p, q][1] == 0, 
      F[p, q][1] == 0}, {p, -.1}, {q, .18}, Method -> "Secant"]


    Out[]= {p -> -0.118505, q -> 0.187442}

Using these roots we find the solution to the original system of equations


bcs = {h[eps] == 1, a[eps] == 0, f'[eps] == 0, 
   f[eps] == -0.1185047596657856`, G[eps] == 0.18744201629754556`};
bcs1 = {G[1] == 0, f[1] == 0};
sol = NDSolve[{ode, bcs}, {a, h, f, G}, {r, eps, 1}, 
   Method -> {"StiffnessSwitching", "NonstiffTest" -> False}];



{Plot[Evaluate[Re[{h[r], h'[r]}] /. sol], {r, eps, 1}, 
  PlotLegends -> {"H", "H'"}, PlotRange -> All], 
 Plot[Evaluate[Re[{a[r], a'[r]}] /. sol], {r, eps, 1}, 
  PlotLegends -> {"A", "A'"}, PlotRange -> All], 
 Plot[Evaluate[Re[{G[r], G'[r]}] /. sol], {r, eps, 1}, 
  PlotLegends -> {"G", "G'"}, PlotRange -> All], 
 Plot[Evaluate[Re[{f[r], f'[r]}] /. sol], {r, eps, 1}, 
  PlotLegends -> {"F", "F'"}, PlotRange -> All]}

fig4 If we want to additionally find L for which h[1]=0, we should move gradually changing the parameters L,p,q.



Comments

Post a Comment

Popular posts from this blog

plotting - How to draw lines between specified dots on ListPlot?

I would like to create a plot where I have unconnected dots and some connected. So far, I have figured out how to draw the dots. My code is the following: ListPlot[{{1, 1}, {2, 2}, {3, 3}, {4, 4}, {1, 4}, {2, 5}, {3, 6}, {4, 7}, {1, 7}, {2, 8}, {3, 9}, {4, 10}, {1, 10}, {2, 11}, {3, 12}, {4,13}, {2.5, 7}}, Ticks -> {{1, 2, 3, 4}, None}, AxesStyle -> Thin, TicksStyle -> Directive[Black, Bold, 12], Mesh -> Full] I have thought using ListLinePlot command, but I don't know how to specify to the command to draw only selected lines between the dots. Do have any suggestions/hints on how to do that? Thank you. Answer One possibility would be to use Epilog with Line : ListPlot[ {{1, 1}, {2, 2}, {3, 3}, {4, 4}, {1, 4}, {2, 5}, {3, 6}, {4, 7}, {1, 7}, {2, 8}, {3, 9}, {4, 10}, {1, 10}, {2, 11}, {3, 12}, {4, 13}, {2.5, 7}}, Ticks -> {{1, 2, 3, 4}, None}, AxesStyle -> Thin, TicksStyle -> Directive[Black, Bold, 12], Mesh -> Full, Epilog -> { Line[ ...
Post a Comment
Read more

dynamic - How can I make a clickable ArrayPlot that returns input?

I would like to create a dynamic ArrayPlot so that the rectangles, when clicked, provide the input. Can I use ArrayPlot for this? Or is there something else I should have to use? Answer ArrayPlot is much more than just a simple array like Grid : it represents a ranged 2D dataset, and its visualization can be finetuned by options like DataReversed and DataRange . These features make it quite complicated to reproduce the same layout and order with Grid . Here I offer AnnotatedArrayPlot which comes in handy when your dataset is more than just a flat 2D array. The dynamic interface allows highlighting individual cells and possibly interacting with them. AnnotatedArrayPlot works the same way as ArrayPlot and accepts the same options plus Enabled , HighlightCoordinates , HighlightStyle and HighlightElementFunction . data = {{Missing["HasSomeMoreData"], GrayLevel[ 1], {RGBColor[0, 1, 1], RGBColor[0, 0, 1], GrayLevel[1]}, RGBColor[0, 1, 0]}, {GrayLevel[0], GrayLevel...
Post a Comment
Read more

Is there a way to do conditional matrix loop using 'continue'

I have the following: n = 3; m = 5; ww = RandomReal[{0, 0.1}, {n, n}]; uu = RandomReal[{0, 1}, {m, n}]; pp = RandomReal[{0, 1}, {n, n}]; ss = RandomInteger[{0, 5}, {m, n}]; Grid[{{"ww", "uu", "pp", "ss"}, {ww // TableForm, uu // TableForm, pp // TableForm, ss // TableForm}}, Spacings -> {5, 2}, Dividers -> All] where I would like to look at every element of matrix ss and produce a matrix tt , with zeroes at the locations in ss which have zeroes, and in all other positions do the following: tt = (-1/Subscript[ww, m]) Log[(1 - uu)/(Subscript[pp, m - 1])], where Subscript[ww, m] is the value at index of ww matrix and where Subscript[pp, m - 1] is the value at index-1 of pp matrix. So for example if the first value ever read from matrix ss happens to be 2, then value taken from matrix ww would be from the row 2, but from pp would be from row 1. Also how to tell difference between a 0 as a valid value from within the matrix elemen...
Post a Comment
Read more