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

front end - keyboard shortcut to invoke Insert new matrix

I frequently need to type in some matrices, and the menu command Insert > Table/Matrix > New... allows matrices with lines drawn between columns and rows, which is very helpful. I would like to make a keyboard shortcut for it, but cannot find the relevant frontend token command (4209405) for it. Since the FullForm[] and InputForm[] of matrices with lines drawn between rows and columns is the same as those without lines, it's hard to do this via 3rd party system-wide text expanders (e.g. autohotkey or atext on mac). How does one assign a keyboard shortcut for the menu item Insert > Table/Matrix > New... , preferably using only mathematica? Thanks! Answer In the MenuSetup.tr (for linux located in the $InstallationDirectory/SystemFiles/FrontEnd/TextResources/X/ directory), I changed the line MenuItem["&New...", "CreateGridBoxDialog"] to read MenuItem["&New...", "CreateGridBoxDialog", MenuKey["m", Modifiers-...
Post a Comment
Read more

How to thread a list

I have data in format data = {{a1, a2}, {b1, b2}, {c1, c2}, {d1, d2}} Tableform: I want to thread it to : tdata = {{{a1, b1}, {a2, b2}}, {{a1, c1}, {a2, c2}}, {{a1, d1}, {a2, d2}}} Tableform: And I would like to do better then pseudofunction[n_] := Transpose[{data2[[1]], data2[[n]]}]; SetAttributes[pseudofunction, Listable]; Range[2, 4] // pseudofunction Here is my benchmark data, where data3 is normal sample of real data. data3 = Drop[ExcelWorkBook[[Column1 ;; Column4]], None, 1]; data2 = {a #, b #, c #, d #} & /@ Range[1, 10^5]; data = RandomReal[{0, 1}, {10^6, 4}]; Here is my benchmark code kptnw[list_] := Transpose[{Table[First@#, {Length@# - 1}], Rest@#}, {3, 1, 2}] &@list kptnw2[list_] := Transpose[{ConstantArray[First@#, Length@# - 1], Rest@#}, {3, 1, 2}] &@list OleksandrR[list_] := Flatten[Outer[List, List@First[list], Rest[list], 1], {{2}, {1, 4}}] paradox2[list_] := Partition[Riffle[list[[1]], #], 2] & /@ Drop[list, 1] RM[list_] := FoldList[Transpose[{First@li...
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