Topic: Differential Equations Boundary-Value Problems
Subject: Numerical Methods
Tool: Scilab
- The Shooting Method
The shooting method is a numerical technique of solving
boundary-value problems (BVPs) for second-order ordinary differential equations
(ODE) of the form shown in Equation (1).
$$y''(x) = f(x, y, y') , \quad a \leq x \leq b \tag{1}$$
with boundary conditions
$$y(a) = \alpha , \quad y(b) = \beta \tag{2}$$
A sample ODE, along with the boundary
conditions provided in Equations (3) and (4), will be solved to demonstrate the
method in this article.
$$\frac{d^2 y}{dx^2} + e^{-x y} + \sin\left(\frac{dy}{dx}\right) = 0 , \quad 1 \leq x \leq 2 \tag{3}$$
$$y(1) = 0 , \quad y(2) = 0 \tag{4}$$
The idea behind the shooting method is to convert the BVP to
the initial-value problem (IVP) as in Equation (5) and (6). The given 2nd-order
ODE requires to be converted to the system of ODEs.
$$y'' = f(x, y, y') \Rightarrow \text{Sys of ODE} \left\{ \begin{aligned} y_1 &= y \\ y_2 &= y' \\ y_1^{'} &= y_2 \\ y_2^{'} &= f(x, y_1, y_2) \end{aligned} \right\} \tag{5}$$
$$y_1(a) = \alpha , \quad y_2(b) = s \tag{6}$$
Solving the converted IVP requires providing a guess for the value of \(s\) and solve the system numerically using Runge-Kutta
method to get the value of the solution, \(y_{\beta} = y(b)\), at \(x = b\). If \(y_{\beta} \neq \beta\), adjust the guess, \(s\), so that \(y_{\beta} - \beta = 0\). For the given two
guesses of \(s\), compute the best guess to make \(y_{\beta} - \beta = 0\) by employing
the root-finding algorithm which is the Bisection method for this article.
For our given problem in Equation (3) and
(4), two guesses value for \(s\) will be tried for \(y_2(a)\).
$$s_1 = 0.1 , \quad s_2 = \frac{\pi}{2} \tag{7}$$
Figure 1 is the resulting solution for
the IVP using \(s_1\) and \(s_2\).
Figure 1. Solution to ODE IVP with \(y(a) = \alpha\) , \(s_1\) , and \(s_2\).
As shown in Figure 1, the two guesses, \(s_1\) and \(s_2\), do not
meet the \(y(b) = \beta\) or \(y(2) = 0\) boundary conditions. Finding the value of \(s\) that makes \(y_{\beta} - \beta = 0\) requires defining a Residual function in Equation
8.
$$\phi(s) = y(b; s) - \beta \tag{8}$$
Residual function measures how far the solution with
guesses is from satisfying the boundary condition at \(x = b\). For the given sample
problem, \(\phi(s) = 1 \times 10^{-6}\) will be utilized in the root-finding
algorithm, the Bisection method. Figure 2 is the visualization of the solution
for sample ODE BVP.
Figure 2. Solution to ODE BVP using the
Shooting Method.
The solution \(y(x)\) satisfies the boundary
conditions \(y(a) = \alpha\) and \(y(b) = \beta\). Implementation of the Runge-Kutta
method and the Bisection method is used in our Scilab script. The two guesses, \(s_1\) and \(s_2\), and the Residual function \(\phi(s)\) are used as parameters in the
script.
- Scilab Script for Figure 2
// Copyright (C) 2025 - I.S.Comia // Date of creation: Jul 12, 2025 // Shooting Method in Solving 2nd-Order ODE BVP clear;clc;clf; // Bundary Value Problem // y'' = -exp(-x.*y)-sin(y'), y(1) = 0, y(2) = 0 // y'' = f(x,y,y'), y(a) = alpha, y(b) = beta // Subroutine // y'' = -exp(-x.*y)-sin(y') in Scilab function function dydx=ode_system(x, y) // y(1) = y, y(2) = y' dydx = [y(2); -exp(x.*y(1))-sin(y(2))]; endfunction // solving the ODE using IVP with guessed initial slope // using Runge-Kutta Method function yb=solve_ivp(s, a, b, h) // Initial value for y and y' x = a:h:b; n = length(x); y = zeros(2, n); y(:, 1) = [0; s]; // y(a) = 0, y'(a) = s (Note: s - slope) // this is Runge-Kutta for i = 1:n-1 k1 = h * ode_system(x(i), y(:, i)); k2 = h * ode_system(x(i) + h/2, y(:, i) + k1/2); k3 = h * ode_system(x(i) + h/2, y(:, i) + k2/2); k4 = h * ode_system(x(i) + h, y(:, i) + k3); y(:, i+1) = y(:, i) + (k1 + 2*k2 + 2*k3 + k4)/6; end yb = y(1, $); // return y(b) endfunction // Main Function // Parameters a = 1; b = 2; alpha = 0; beta = 0; h = 0.05; // Two Initial guesses for s = y'(a) or slopes, s @ x = a s1 = 0.1; s2 = %pi/2; tol = 1e-6; // Plot of IVP using s1 and s2 x = a:h:b x0 = a y1_0 = [0;s1] y1 = ode(y1_0,x0,x,ode_system) y2_0 = [0;s2] y2 = ode(y2_0,x0,x,ode_system) plot(x, y1(1, :), "g.-", "linewidth",1); plot(x, y2(1, :), "g.-", "linewidth",1); // Shooting Method with Bisection Method while abs(s2 - s1) > tol smid = (s1 + s2)/2; yb1 = solve_ivp(s1, a, b, h); ybmid = solve_ivp(smid, a, b, h); if (yb1 - beta)*(ybmid - beta) < 0 then s2 = smid; else s1 = smid; end end // Final solution with best guess s_final = (s1 + s2)/2; x = a:h:b; n = length(x); y = zeros(2, n); y(:, 1) = [alpha; s_final]; for i = 1:n-1 k1 = h * ode_system(x(i), y(:, i)); k2 = h * ode_system(x(i) + h/2, y(:, i) + k1/2); k3 = h * ode_system(x(i) + h/2, y(:, i) + k2/2); k4 = h * ode_system(x(i) + h, y(:, i) + k3); y(:, i+1) = y(:, i) + (k1 + 2*k2 + 2*k3 + k4)/6; end // Visualization of solution to ODE BVP fig = gcf() fig.figure_size = [700, 700] plot(x, y(1, :), "b-", "linewidth",3); plot(x, y(1, :), "ro", "linewidth",4.5); title("$\text{Shooting Method Solution to ODE BVP}$","fontsize",4.5); xlabel("$\LARGE x$") ylabel("$\LARGE y(x)$") note = ["$\Large y(x) \; \text{w/ IVPs}$";"$\Large y(x) \; \text{w/ BVP}$";" "] legend(note,2,with_box=%F) xgrid(color("grey"),1,7) // Plot annotation h1 = xstring(1.55,0.725,"https://gani-mech-toolbox.blogspot.com") h1.font_size = 2 eqn1 = "$\frac{d^2 y}{dx^2}+e^{-xy}+\sin\left(\frac{dy}{dx}\right)=0\text{,}\;
y(1)=0\text{,}\;y(2)=0$" h2 = xstring(1.05,-0.4,eqn1) h2.font_size = 3.5 h3 = xstring(1.05,-0.25,"Given : 2nd-Order ODE BVP") h3.font_size = 3 eqn2 = "$y(a) = \alpha, \quad \left.\frac{dy}{dx} \right|_{x=a} = s_1$" h4 = xstring(1.6,-0.125,eqn2) h4.font_size = 3.5 eqn3 = "$y(a) = \alpha, \quad \left.\frac{dy}{dx} \right|_{x=a} = s_2$" h5 = xstring(1.5,0.55,eqn3) h5.font_size = 3.5 eqn4 = "$y(a) = \alpha, \quad \left.\frac{dy}{dx} \right|_{x=a} \; \mid \;
y_\beta - \beta = 0$" h6 = xstring(1.5,0.15,eqn4) h6.font_size = 3.5 ax = gca() ax.data_bounds = [1 -0.4; 2 0.8]
- Scilab Script for Figure 1
// Copyright (C) 2025 - I.S.Comia // Date of creation: Jul 12, 2025 // Solution to 2nd-Order ODE IVP using Scilab's ode() function clear;clc;clf; // Initial Value Problem // y'' = -exp(-x.*y)-sin(y'), y(1) = 0, y'(1) = s1 & s2 // y'' = f(x,y,y'), y(a) = alpha, y'(a) = s1 & s2 // Subroutine // y'' = -exp(-x.*y)-sin(y') in Scilab function function dydx=ode_system(x, y) // y(1) = y, y(2) = y' dydx = [y(2); -exp(x.*y(1))-sin(y(2))]; endfunction // Main Function // Parameters a = 1; b = 2; alpha = 0; h = 0.05; // Two Initial guesses for s = y'(a) or slope s @ x = a s1 = 0.1; s2 = %pi/2; // Plot of IVP using s1 and s2 x = a:h:b x0 = a y1_0 = [0;s1] y1 = ode(y1_0,x0,x,ode_system) y2_0 = [0;s2] y2 = ode(y2_0,x0,x,ode_system) // Visualization of solution to ODE IVP fig = gcf() fig.figure_size = [700, 700] plot(x, y1(1, :), "g.-", "linewidth",1); plot(x, y2(1, :), "m.-", "linewidth",1); note = ["Solution to ODE IVP with ,","$y(a)=\alpha$",",","$\frac{dy}{dx}(a) = s$"] title(note,"fontsize",3.5); xlabel("$\LARGE x$") ylabel("$\LARGE y(x)$") note1 = ["$y(x) \; \text{w/ IVP} \; s_1 = 0.1$","$y(x) \; \text{w/ IVP} \;
s_2 = \frac{\pi}{2}$"] a1 = legend(note1,2,with_box=%F) a1.font_size = 3 xgrid(color("grey"),1,7) // Plot annotations h1 = xstring(1.55,0.725,"https://gani-mech-toolbox.blogspot.com") h1.font_size = 2 eqn1 = "$\frac{d^2 y}{dx^2}+e^{-xy}+\sin\left(\frac{dy}{dx}\right)=0\text{,}\;
y(1)=0\text{,}\;\frac{dy}{dx}(1)=s$" h2 = xstring(1.05,-0.4,eqn1) h2.font_size = 3.5 h3 = xstring(1.05,-0.25,"Given : 2nd-Order ODE IVP") h3.font_size = 3 eqn2 = "$y(a) = \alpha, \quad \left.\frac{dy}{dx} \right|_{x=a} = s_1$" h4 = xstring(1.6,-0.125,eqn2) h4.font_size = 3.5 eqn3 = "$y(a) = \alpha, \quad \left.\frac{dy}{dx} \right|_{x=a} = s_2$" h5 = xstring(1.5,0.55,eqn3) h5.font_size = 3.5 // Plot limits ax = gca() ax.data_bounds = [1 -0.4; 2 0.8]
The shooting method is an effective
approach for solving second-order boundary value problems (BVPs), particularly
in mechanical, thermal, and structural engineering applications. To apply this
method, the BVP is converted into an initial value problem (IVP) by estimating
the necessary initial conditions. The IVP is then solved, and the results are
compared with the boundary conditions at the opposite end. This process is
repeated, adjusting the initial guesses as needed, until the solution satisfies
the boundary conditions. While this approach involves some trial and error,
linear second-order problems typically require no more than two iterations to
converge.
Feel free to comment for inquiry,
clarification, correction or suggestion for improvement. Drop your email to make
a request to the author.
References
- Gerald, C. and Wheatly, P., Applied Numerical Analysis. Addison-Wesley Publishing Company, Inc. 1994.
- ChatGPT. AI Tools. Accessed on July 12, 2025.
- Google AI Overview. AI Tools. Accessed on July 12, 2025.
