A Note on the Shooting Method

Allan Hayes

Mathematica in Education
Vol. 3 No.4
Copyright 1994
TELOS/Springer-Verlag

In Volume 2 Number 4, Robert Spall gave an interesting application of Mathematica to fluid dynamics. The mathematical part was concerned with the solution of a  nonlinear 2-point boundary by the shooting method. The present note takes the reduction of programming effort a bit further by using FindRoot instead of a While loop. Some printout about progress is lost but the payoff is simpler code and better display of the mathematical ideas.

The system to be solved is .

{f[r](f''[r] - 1/r f'[r] + 
		2 V0^2/U Exp[-2/U f[r]] (1.0 - Exp[-2/U f[r]])) +
		V0^2 r^2 (Exp[-4/U f[r] ] - Exp[-2/U f[r]]) == 0,
		f[b1] == 0.5(a1^2 U),
		f[b2] == d
		};

We fix the parameters  other than V0.

U = 1.0;
a1 = 4.0; b1 = 4.1;
a2 = .1^10;  b2 = .1^10;
d = 0.5(a2^2 U);

The system has two boundary points b1 and b2 with conditions  f[b1] == 0.5(a1^2 U) and f[b2] == d . NDSolve cannot solve such a system.  The crucial idea is to replace the boundary condition f[b2] == d at b2 with the variable condition f'[b1] == s at b1, so that we have a one boundary point system. For each numerical value of s we have a one  NDSolve can solve this system. It remains for us to find a value s0 of s  such that corresponding solution f0 of the modified system has f0[b2] = d , making it a solution of the original equation. One more detailed technical point, since we  want to use f' we introduce a symbol for it so that it can be calculated by NDSolve as the equation is being solved. Now we start the solution. Here is the new system with its two parameters displayed and accessible.

sys[s_,V0_] := 
	{f[r](f1'[r] - 1/r f1[r] + 
		2 V0^2/U Exp[-2/U f[r]] (1.0 - Exp[-2/U f[r]])) +
		V0^2 r^2 (Exp[-4/U f[r] ] - Exp[-2/U f[r]]) == 0,
		f'[r] == f1[r],
		f[b1] == 0.5(a1^2 U),
		f1[b1] == s
	};

The following function returns  the pair {f, f1} for different  s and V0.

g[s_,V0_] := {f,f1} /. NDSolve[sys[s,V0],{f,f1},{r,b2,b1}][[1]]

For each value of V0 the following function finds a value of s that makes f a solution of the original equation and remembers pair {f,f1} as its value. Notice that we have to give FindRoot two initial guesses b1 U and 1.05 b1 U (please refer to the original article for the reason for choosing these).

sol[V0_] :=
	(	FindRoot[
			(sol[V0] = g[s,V0])[[1]][b2] == d, {s, b1 U, 1.05 b1 U},
			AccuracyGoal -> 10^-8
		];
		sol[V0]		
	)

We can plot graphs simultaneously.

multiplot[V0___] :=
	Plot[ Evaluate[(1/r sol[#][[2]][r])&/@{V0}], {r,b2,b1},
		Frame -> True,
		FrameLabel -> {"r", "u"},
		PlotRange -> {-1.25, 1.25},
		Axes -> False,
		Prolog -> {{RGBColor[1,0,0],Line[{{0,1},{4.1,1}}]}}
	]

For example:

multiplot[0.0,0.49, 1.0,1.2, 1.5]

Or singly

plot[V0_] :=
	Plot[ Evaluate[1/r sol[V0][[2]][r]], {r,b2,b1},
		Frame -> True,
		FrameLabel -> {"r", "u"},
		PlotRange -> {-1.25, 1.25},
		Axes -> False,
		PlotLabel -> StringForm["V0 = ``",V0],
		Prolog -> {{RGBColor[1,0,0],Line[{{0,1},{4.1,1}}]}}
	]

For example:

plot[1.2]

And make animations.

For example:

Do[plot[V0],{V0, 0, 1.5, .1}]