Shooting Method with extra Unknown Condition

Question

I've been trying to solve a system of 3 coupled 2nd order ODEs, for a real variable $x$, $0geq xleq infty$. The equations are the following:
begin{align}
&x^{2},h''(x) - x,h'(x) + x^{2},g^{2}(x)left[1-h(x)right] = 0,,
&x^{2},f''(x) + x,f'(x) - lambda, x^{2},f(x)left[f^{2}(x) + g^{2}(x) - 2right] = 0,,
&x^{2},g''(x) + x,g'(x) - frac{1}{2},g(x)left[1-h(x)right]^{2} - lambda, x^{2}g(x)left[f^{2}(x) + g^{2}(x) - 2right] = 0,.
end{align}
In addition, the BCs are (where my problem starts):
$$h(0)=0=g(0),, quad f(0)=Omega$$
and
$$h(xtoinfty)=f(xtoinfty)=g(xtoinfty)=1,.$$
First of all, I decided to solve to some finite $x$ such as $x_{max}$ and then try to increase this domain. Then, my problem consists of how I should "tell" NDSolve that I don't know the value $Omega$ is going to have. I know that, somehow, the numerical solution must find an appropriate value for $Omega$ that agrees with the whole solution. But I cannot understand how I can do this. My starting code is the following:
lambda = 0.5; 
eps = 0.001;
xmax = 5;
eq1=x^2*h''[x] - x*h'[x] + x^2*(g[x]^2) (1 - h[x]);
eq2= x^2*f''[x] + x*f'[x] - lambda*x^2*f[x] ((f[x]^2) + (g[x]^2) - 2);
eq3= x^2*g''[x] + x*g'[x] -  1/2*g[x] (1 - h[x])^2 - lambda*x^2*g[x] ((f[x]^2) + (g[x]^2) - 
2);
    
sols=First[NDSolve[{eq1==0,eq2==0, eq3== 0,h[eps] == 0,f[eps] == Omega, g[eps]==0}, {f[x], 
g[x], h[x]}, {x, eps, xmax},Method -> {"Shooting","StartingInitialConditions" -> {h[eps] == 
0,f[eps] == Omega, g[eps] == 0}}, WorkingPrecision -> 5]];

As you can see, my code is incomplete. The shooting method would need 6 initial conditions for the (converted) IVP -> 3 from the BC at $x=0$ and the shooting for the 3 first-order derivatives. However, since I don't know (a priori) the value of $Omega$, I'm stuck :(
Ps.: From my problem, I know I can put by hand that all the first-order derivatives go to zero when $xtoinfty$... But this would add too many conditions for Mathematica, right?
Could you, please, provide any advice on how I can tackle the problem?

bbgodfrey · Accepted Answer

With NDSolve and Method -> "Shooting", I was unable to integrate past x = 2.93. With Method -> "FiniteElement", I was able to integrate as far as x = 7.5, but omega could only be approximated. In general, the problem appeared to be inadequate numerical accuracy near x = 0. The following worked much better. First, obtain approximate symbolic solutions near x = 0 sh0 = DSolveValue[{eq1 == 0, h[0] == 0} /. {f[x]^2 -> omega^2, g[x]^2 -> 0}, h[x], x] /. C[1] -> ch (* (ch x^2)/2 *) sf0 = Simplify[DSolveValue[{eq2 == 0, f[0] == omega} /. {f[x]^2 -> omega^2, g[x]^2 -> 0}, f[x], x], omega^2 < 2] (* omega BesselJ[0, (Sqrt[2 - omega^2] x)/Sqrt[2]] *) sg0 = Simplify[DSolveValue[{eq3 == 0, g[0] == 0} /. {f[x]^2 -> omega^2, g[x]^2 -> 0, h[x] -> 0}, g[x], x], omega^2 < 2] /. C[1] -> cg (* cg (BesselJ[1/Sqrt[2], (Sqrt[2 - omega^2] x)/Sqrt[2]] - (BesselJ[1/Sqrt[2], 0] BesselY[1/Sqrt[2], (Sqrt[2 - omega^2] x)/Sqrt[2]])/BesselY[1/Sqrt[2], 0]) *) Then, obtain the three parameters {omega, cg, ch} by what might be called do-it-yourself shooting, with the symbolic solutions near x = 0 as initial conditions. xmax = 10; eps = .1; sp = ParametricNDSolveValue[{eq1 == 0, eq2 == 0, eq3 == 0, {f[x] == sf0, g[x] == sg0, h[x] == sh0, f'[x] == D[sf0, x], g'[x] == D[sg0, x], h'[x] == D[sh0, x]} /. x -> eps}, {f[xmax], g[xmax], h[xmax]}, {x, eps, xmax}, {omega, cg, ch}, PrecisionGoal -> 10, AccuracyGoal -> 10]; FindRoot[sp[omega, cg, ch] - 1, {{omega, 1.35617}, {cg, 1.3415}, {ch, 0.325786}}, Evaluated -> False] (* {omega -> 1.35617, cg -> 1.3415, ch -> 0.325786} *) Finally, compute and solve the equations with the parameters just determined. NDSolveValue[{eq1 == 0, eq2 == 0, eq3 == 0, {f[x] == sf0, g[x] == sg0, h[x] == sh0, f'[x] == D[sf0, x], g'[x] == D[sg0, x], h'[x] == D[sh0, x]} /. x -> eps} /. %, {f[x], g[x], h[x]}, {x, eps, xmax}]; Plot[%, {x, eps, xmax}, ImageSize -> Large, AxesLabel -> {x, "f,g,h"}, LabelStyle -> {15, Bold, Black}] The initial guesses for FindRoot were obtained by integrating the equations for xmax = 3 and using the result as guesses for xmax = 4, etc. Note that even xmax = 10 is not in the asymptotic domain of the equations. Increasing xmax to, say 20 undoubtedly would require higher WorkingPrecision and an automated process for gradually increasing xmax, both of which are feasible.

Shooting Method with extra Unknown Condition

One Answer

Add your own answers!

Ask a Question