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

equation solving - Invert and fit implicitly defined curve

I need to fit an implicitly defined curve. I thought I could get some data out of Solve , and then using FindFit . Therefore, I would like to find the relation the parametric curve defined by $F(x,y)=0$: Solve[-(1/2) + 1/2 (0.41202 BesselK[0, 0.1 Sqrt[x^2 + y^2]] + (0.101483 x BesselK[1, 0.1 Sqrt[x^2 + y^2]])/Sqrt[x^2 + y^2]) == 0, y] But I can't get an output: Solve was unable to solve the system with inexact coefficients or the system obtained by direct rationalization of inexact numbers present in the system. Since many of the methods used by Solve require exact input, providing Solve with an exact version of the system may help. >> Edit: In particular, I would like to fit the data coming from the curve with the expression of another curve, and not with a function $f(x)$. In particular, since this clearly looks like a cardioid , I would like it to fit to something like it. What other strategies could I try?
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