von Mises Stress Failure Theory in a Spinning Solid Disk

Topic: von Mises Yield Criterion and Rotating Disk

Subject: Machine Design

Tool: Scilab

By: Gani Comia, September 2025

• Spinning Solid Disks

Rotating circular machine elements can be modeled as rotating disks to evaluate stresses. Similar to the theory of thick-walled cylinders, both tangential and radial stresses are present, but in this case, they result from inertial forces within the ring or disk. The tangential and radial stresses are governed by the following conditions.

1. The outside radius of the disk is significantly larger than its thickness, \(r_o \geq 10 \, t\).
2. The ring or disk has a uniform thickness.
3. The stresses are evenly distributed throughout the thickness.

Reference article and derivation of the stresses in the rotating rings and disks can be found in the website https://roymech.org/Useful_Tables/Mechanics/Rotating_cylinders.html.

• Tangential Stress Equation for Solid Disk

$$\sigma_t = \left(\frac{\rho \, \omega^2}{8}\right) \left[ (3 + \nu) \,  r_o^2 - (1 - 3 \, \nu) \, r^2 \right] \tag{1}$$

Where:

\(\sigma_t\) – tangential stress, \(Pa\)

\(\rho\) - density, \(kg/m^2\)

\(\omega\) – angular velocity, \(rad/s\)

\(\nu\) – Poisson’s Ratio

\(r_o\) – outside radius of cylinder, \(m\)

\(r\) – radius at the point of calculation, \(m\)

• Radial Stress Equation for Solid Disk

$$\sigma_r = \left( \frac{\rho \, \omega^2}{8} \right) (3 + \nu) (r_o^2 - r^2) \tag{2}$$

Where:

\(\sigma_r\) – radial stress, \(Pa\)

• von Mises Yield Criterion (Maximum Distortion Energy Criterion)

In engineering, the von Mises yield criterion states that yielding begins when the von Mises stress \(\sigma_v\)​, derived from the Cauchy stress tensor, equals the yield strength \(\sigma_y\)​.

The von Mises Equation for a general state of stress is

$$\sigma_{\nu} = \sqrt{\frac{1}{2} \left[ (\sigma_{11} - \sigma_{22})^2 + (\sigma_{22} - \sigma_{33})^2 + (\sigma_{33} - \sigma_{11})^2 \right] + 3 \, (\sigma_{12}^2 + \sigma_{23}^2 + \sigma_{31}^2)} \tag{3}$$

For the boundary conditions,

$$\sigma_{12} = \sigma_{23} = \sigma_{31} = 0 \tag{4}$$

the von Mises defines the equations under principal stress

$$\sigma_{\nu} = \sqrt{\frac{1}{2} \left[ (\sigma_1 - \sigma_2)^2 + (\sigma_2 - \sigma_3)^2 + (\sigma_3 - \sigma_1)^2 \right]} \tag{5}$$

For the 2D multi-axial stress, where \(\sigma_3 = 0\) as boundary condition under principal plane stress, the equivalent von Mises stress, \(\sigma_{\nu}\), becomes

$$\sigma_{\nu} = \sqrt{\sigma_1^2 + \sigma_2^2 - \sigma_1 \, \sigma_2} \tag{6}$$

• von Mises Yield Criterion and Factor of safety

The use of the von Mises criterion as a yield criterion is only exactly applicable when the material properties are isotropic, there is no considerable dynamic or unpredictable loads, and the ratio of the shear yield strength to the tensile yield strength has the value as in Equation (7).

$$\frac{\sigma_s}{\sigma_y} = \frac{1}{\sqrt{3}} \approx 0.577 \tag{7}$$

In practice it is necessary to use engineering judgement and use the above ratio.

The factor of safety based on yield strength is given by Equation (8).

$$N_y = \frac{\sigma_y}{\sigma_{\nu}} - 1 \tag{8}$$

Critical components suggest a factor of safety of \(\text{2.5}  \leq  N_y  \leq  \text{3.0} \).

• Application

A solid disk (Figure 1) is analyzed at different rotational speeds to determine safe operating limits. The design factor indicates the maximum speed at which the disk can rotate without failure. Using the von Mises Yield Criterion, the equivalent stress is evaluated from the combined tangential and radial stresses. The design factor, as given by Equation 8, defines the likelihood of structural failure for the specified parameters.

Assumed Parameters:

Material:    Alloy Steel

Density:    7,850 \(kg/m^3\)

Yield Strength:    1,250 \(MPa\)

Poisson’s Ratio:    0.29

Radius:    0.5 \(m\)

Rotational Speed:    2500, 5000, 7500, and 10000 \(rpm\)

Figure 1. Rotating Solid Disk (ChatGPT Image Sep 10, 2025).

Figure 2 presents the equivalent combined stresses according to the von Mises yield criterion. The stress distribution varies radially across the disk, reaching its maximum at the center where \(r = 0\).

Figure 2. von Mises Stress from the Point of Analysis along the Disk Radius.

A contour plot, shown in Figure 3, provides an alternative visualization of the stress distribution for a given rotational speed.

Figure 3. von Mises Stress Distribution

Figure 4 shows the design factor at various rotational speeds. The minimum factor of safety is indicated for each case. Rotational speed of \(7500  \, \text{rpm} \) and above do not satisfy the criteria for critical components.

Figure 4. Design Factor using von Mises Stress for Spinning Disk.

Refer to the following Scilab scripts for figure reproduction.

Scilab Script for Figure 2 and 4

// Copyright (C) 2025 - Gani Comia
// Date of creation: 12 Sep 2025
// Script: von Mises Yield Criterion in a Spinning Disk
clear;clc;

// Primary parameters
rho = 7850;                             // kg/m^3, density of steel
sigmaY = 1250;                          // MPa, yield strength, alloy steel
rpm = [2500 5000 7500 10000];           // rpm, rotational speed      
radps = rpm.*(2*%pi/60);                // rad/s, maximum rotational speed
nu = 0.29;                              // Poisson's ration of carbon steel
r_o = 0.5;                              // m, outer radius of cylinder

// Secondary parameters
R = linspace(0,r_o,100)
n = length(radps)

// von Mises stress distribution for a rotational speed
for i = 1:n
    // tangential stress, Pa
    sigmaT(i,:) = ((rho.*radps(i).^2)/8).*((3+nu).*(r_o).^2 - (1-3*nu).*(R.^2));
    // radial stress, Pa
    sigmaR(i,:) = ((rho.*radps(i).^2)/8).*((3+nu).*((r_o).^2 - (R.^2)));
    // von Mises stress, Pa
    sigmaV(i,:) = sqrt(sigmaT(i,:).^2 + sigmaR(i,:).^2 - (sigmaT(i,:).*sigmaR(i,:)))
    // von Mises stress, MPa
    sigmaV(i,:) = sigmaV(i,:)/1e6
end

// Plot of von Mises stress distribution for a given rotation speed
clf;
scf(0)
f = gcf()
f.figure_size = [800,800]
subplot(2,2,1)
    plot(R,sigmaV(1,:),"g-","linewidth",4.5);
    title("von Mises Stress at 2,500 rpm")
    xlabel("Disk Radius, m","fontsize",3)
    ylabel("Stress, MPa","fontsize",3)
    xgrid(color("grey"),1,7)
    mprintf("Max von Mises #1 : %0.0f\n", max(sigmaV(1,:)))
    legend("$\LARGE \sigma_{Vmax} = \text{55 MPa}$",with_box=%f)
subplot(2,2,2)
    plot(R,sigmaV(2,:),"m-","linewidth",4.5);
    title("von Mises Stress at 5,000 rpm")
    xlabel("Disk Radius, m","fontsize",3)
    ylabel("Stress, MPa","fontsize",3)
    xgrid(color("grey"),1,7)
    mprintf("Max von Mises #3 : %0.0f\n", max(sigmaV(2,:)))
    legend("$\LARGE \sigma_{Vmax} = \text{221 MPa}$",with_box=%f)
subplot(2,2,3)
    plot(R,sigmaV(3,:),"b-","linewidth",4.5);
    title("von Mises Stress at 7,500 rpm")
    xlabel("Disk Radius, m","fontsize",3)
    ylabel("Stress, MPa","fontsize",3)
    xgrid(color("grey"),1,7)
    mprintf("Max von Mises #3 : %0.0f\n", max(sigmaV(3,:)))
    legend("$\LARGE \sigma_{Vmax} = \text{498 MPa}$",with_box=%f)
subplot(2,2,4)
    plot(R,sigmaV(4,:),"r-","linewidth",4.5);
    title("von Mises Stress at 10,000 rpm")
    xlabel("Disk Radius, m","fontsize",3)
    ylabel("Stress, MPa","fontsize",3)
    xgrid(color("grey"),1,7)
    mprintf("Max von Mises #4 : %0.0f\n", max(sigmaV(4,:)))
    legend("$\LARGE \sigma_{Vmax} = \text{885 MPa}$",with_box=%f)
    a1 = xstring(0.25,860,"https://gani-mech-toolbox.blogspot.com")
    a1.font_size = 2
    
// factor-of-safety calculation and plot
for j = 1:n
    Ny(j,:) = (sigmaY./sigmaV(j,:)) - 1
end

scf(1)
g = gcf()
g.figure_size = [700,700]
plot("ln",R,Ny(1,:),"g-","linewidth",4.5)
plot("ln",R,Ny(2,:),"m-","linewidth",4.5)
plot("ln",R,Ny(3,:),"b-","linewidth",4.5)
plot("ln",R,Ny(4,:),"r-","linewidth",4.5)
title("Design Factor using von Mises Stress Criterion for Spinning Disk","fontsize",3.5)
xlabel("Disk Radius, m","fontsize",3.5)
ylabel("Design Factor, Ny","fontsize",3.5)
xgrid(color("grey"),1,7)

// minimum design factor for a given rotation annotation
format(5)
a2 = xstring(0.005,Ny(1,1),["N (min) = ",string(Ny(1,1))," at 2,500 rpm"])
a2.font_size = 3
a3 = xstring(0.005,Ny(2,1),["N (min) = ",string(Ny(2,1))," at 5,000 rpm"])
a3.font_size = 3
a4 = xstring(0.005,Ny(3,1),["N (min) = ",string(Ny(3,1))," at 7,500 rpm"])
a4.font_size = 3
a5 = xstring(0.005,Ny(4,1)-0.2,["N (min) = ",string(Ny(4,1))," at 10,000 rpm"])
a5.font_size = 3

// design parameter annotation
a6 = "$N_y = \frac{\sigma_y}{\sigma_v} - 1$"
a7 = "$\sigma_v = \sqrt{\sigma_1^2 + \sigma_2^2 - \sigma_1 \sigma_2}$"
a8 = "$\frac{\sigma_s}{\sigma_y} = \frac{1}{\sqrt{3}} \; , \; \text{assumption}$"
a9 = xstring(0.0625,10,[a6;a7;a8])
a9.font_size = 4
a11 = xstring(0.05,8,"https://gani-mech-toolbox.blogspot.com")
a11.font_size = 2.5

ax = gca()
ax.data_bounds = [0.005 0; 0.5 25];

Scilab Script for Figure 3

// Copyright (C) 2025 - Gani Comia
// Date of creation: 5 Sep 2025
// Script: von Mises Stress Distribution
clear;clc;

// Primary parameters
rho = 7850;                             // kg/m^3, density of steel
rpm = [2500 5000 7500 10000];           // rpm, rotation speed
radps = max(rpm.*(2*%pi/60));           // rad/s, maximum rotational speed
nu = 0.29;                              // Poisson's ration of carbon steel
r_o = 0.5;                              // m, outer radius of cylinder

// Create meshgrid
n = 250;
x = linspace(-r_o, r_o, n);
y = linspace(-r_o, r_o, n);
[X, Y] = ndgrid(x, y);

// Compute radial distance from center
R = sqrt(X.^2 + Y.^2);

// Define stress distribution at maximum rotational speed
sigmaT = zeros(n, n);
sigmaR = zeros(n, n);
sigmaV = zeros(n, n);

for i = 1:n
    for j = 1:n
        if R(i,j) <= r_o then
            // tangential stress
            sigmaT(i,j) = ((rho.*radps.^2)/8).*((3+nu).*(r_o).^2 - (1-3*nu).*(R(i,j)).^2);
            // radial stress
            sigmaR(i,j) = ((rho.*radps.^2)/8).*((3+nu).*((r_o).^2 - (R(i,j)).^2));
            // von Mises stress
            sigmaV(i,j) = sqrt(sigmaT(i,j).^2 + sigmaR(i,j).^2 - (sigmaT(i,j).*sigmaR(i,j))); 
        else
            sigmaT(i,j) = %nan;     // Outside the circle
            sigmaR(i,j) = %nan;
            sigmaV(i,j) = %nan;
        end
    end
end

// Stress distribution using contour plot
clf;
f = gcf()
f.figure_size = [700,700]

sigmaV = sigmaV./(1e6)
disp(max(sigmaV)); disp(min(sigmaV));

contourf(x, y, sigmaV, 60);
title("von Mises Stress Distribution in Circular Domain at 10K rpm","fontsize",3.5);
xlabel("x-coordinates","fontsize",3)
ylabel("y-coordinates","fontsize",3)
xgrid(color("grey"),1,7)
isoview on;

f.color_map = jet(64)
colorbar(750,900)

ax = gca()
ax.data_bounds = [-0.6 -0.6; 0.6 0.6]
ax.auto_ticks = ["on" "on"]

Feel free to comment for inquiry, clarification, correction or suggestion for improvement. Drop your email to make a request to the author.

Disclaimer: The formulas and calculations presented are for technical reference only. Users must verify the accuracy and ensure compliance with applicable engineering standards, codes, and safety requirements before practical application.

References

1. “Rotating Disks and Cylinders”. RoyMech.org. 2020. https://roymech.org/Useful_Tables/Mechanics/Rotating_cylinders.html .
2. “von Mises yield criterion”. Wikipedia The Free Encyclopedia. 18 September 2024. https://en.wikipedia.org/wiki/Von_Mises_yield_criterion .
3. J.E. Shigley and C.R. Mischke. Mechanical Engineering Design. 5th Ed. New York, New York 10020. McGraw-Hill Publishing Company. 1989. 

Factor of Safety for Basic Design of Machine Elements under Simple Stress

Topic: Design Factor and Simple Stresses

Subject: Machine Design

Tool: Scilab & QCAD

By: Gani Comia, August 2025

• Types of Design

1. Rational design – regarded as purely mathematical and founded on the principles of mechanics.
2. Empirical design – a design approach that follows established practices, relying primarily on past performance rather than theoretical justification.
3. Industrial design – a term generally used for designing with emphasis on both appearance and function.

Design that relies on the factor of safety, such as in cases of simple static compression or tensile loads, is classified as rational design. On the other hand, industrial design addresses aspects like final appearance and additional features of machine elements.

• Simple or Direct Stress

A fundamental challenge for engineers is selecting the appropriate material and applying it correctly and proportionately, ensuring that a structure or machine element performs its intended function efficiently. To achieve this, it is crucial to determine the material’s strength, stiffness, and other properties.

The unit strength of a material is generally defined in terms of stress. Stress is denoted symbolically as:

$$\sigma = \frac{F}{A} \tag{1}$$

Where:

\(\sigma\) – stress or force per unit area

\(F\) – applied compressive or tensile force

\(A\) – cross-sectional area

• Factor of Safety or Design Factor

The Factor of Safety (\(FoS\)) is a value used to divide the strength criterion to a design stress in order to establish a design criterion. In a literal sense, it represents how many times stronger the design is compared to the expected load. However, since the term does not precisely match its literal meaning, it is sometimes referred to as the design factor. The design factor, or factor of safety, establishes the design stress based on ultimate-strength or yield-strength considerations.

$$\sigma_d = \frac{\sigma_u}{N}, \quad \sigma_d = \frac{\sigma_y}{N} \tag{2}$$

Where:
\(\sigma_d\) – design stress

\(\sigma_u\) – ultimate-stress criteria

\(\sigma_y\) – yield-stress criteria
\(N\) – factor of safety or design factor

In this article, since the machine element under consideration is made of ductile material, the following allowable values of design factor should be applied under dead load conditions.

$$\text{Dead Load, } N = \left\{ \begin{aligned} N & \geq 3 \; \; \; \, \text{ based on } \sigma_u \\ N & \geq 1.5 \; \text{ based on } \sigma_y \end{aligned} \right\} \tag{3}$$

The more basic definition of factor of safety for a single load is:

$$\text{Factor of Safety} = \frac{\text{Load that would cause failure}}{\text{Actual load on part}} \tag{4}$$

• Application Example

Figure 1 shows a steel dead weight of \(8 \; \text{tons}\) supported by four cylindrical polyurethane material. Three types of polyurethane materials with ultimate strengths of \(\sigma_u = 10, \; 20, \; \text{and} \; 30 \text{MPa}\) are considered in this example. Each cylindrical support has a diameter of \(4 \; \text{inches}\). Based on this information, the factor of safety will be evaluated for each material type under a dead load ranging from \(1 \; \text{to} \; 12 \; \text{tons}\).

Figure 1. Simply Supported Dead Load.

One of the key concepts in mechanics is the free-body diagram (FBD). A FBD is a sketch of an isolated body that illustrates the external forces acting on it. The reaction forces represent the forces exerted by the body on its supports or adjacent bodies. Figure 2 shows the free-body diagram of a simply supported dead load.

Figure 2. Free-body Diagram for Simply Supported Dead Load of 8 Tons.

Analysis of the given scenario are shown in Figure 3. Three different factors of safety is shown on the figure for each different materials defined by its ultimate strength, \(\sigma_u\), at a particular \(8 \; \text{tons}\) weight or load.

Figure 3. Factor of Safety Calculation for a Simple Stress Analysis.

Visualizing the graph of the calculation provides far more insight than a single computed value for a given scenario. The Scilab script below recreates the plot shown in Figure 3.

• Scilab Script

// Copyright (C) 2025 - Gani Comia
// Date of creation: 24 Aug 2025
// Factor-of-Safety Calculation for Simple Compressive Stress
clear;clc;clf;

// primary parameters
W = linspace(1,12,50);                  // tons, compressive load or weight
Ncyl = 4;                               // pcs, number of cylindrical support
F = 9806.65*(W/Ncyl);                   // N, load
d = 4-0.02;                             // inch, diameter (min)
d = 25.4*d;                             // mm, diameter (min)
A = (%pi/4)*d.^2;                       // mm^2, area

// sample material properties of a brittle material
sigma = [10 20 30];                    // MPa, ultimate or yield strength
Nst =  length(sigma)

// calculation of factor-of-safety for different material type
sigmaD = F./A
for i = 1:Nst
    FoS(i,:) = sigma(i)./sigmaD
end

// visualization
f = gcf();
f.figure_size = [700,700];
plot(W,FoS(1,:),"b-","linewidth",3)
plot(W,FoS(2,:),"c-","linewidth",3)
plot(W,FoS(3,:),"g-","linewidth",3)
title("Factor-of-Safety for Simple Compressive Stress","fontsize",3.5)
xlabel("Weight, W, (tons)","fontsize",3.25)
ylabel("Factor-of-Safety @ Cyl Support, FoS","fontsize",3.25)
note1 = "$\text{Matl#1},\;\sigma_u = 10\;\text{MPa}$"
note2 = "$\text{Matl#2},\;\sigma_u = 20\;\text{MPa}$"
note3 = "$\text{Matl#3},\;\sigma_u = 30\;\text{MPa}$"
leg = legend([note1,note2, note3],with_box=%F)
leg.font_size = 3
note2 = xstring(6.5,23-0.25,"https://gani-mech-toolbox.blogspot.com")
note2.font_size = 2
xgrid(color("grey"),1,7)

ax = gca()
ax.data_bounds = [0 0; 12 30]

// plotting intersecting lines & markers
idW = max(find(W <= 8.0))
W_val = W(idW)
mprintf("W_val: %3.1f, FoS_val: %3.1f \n",W_val,FoS(1,idW))
mprintf("W_val: %3.1f, FoS_val: %3.1f \n",W_val,FoS(2,idW))
mprintf("W_val: %3.1f, FoS_val: %3.1f \n",W_val,FoS(3,idW))

W_line_v = [8 8]; FoS_line_v = [0 FoS(3,idW)];
plot(W_line_v,FoS_line_v,"r--","linewidth",1)

for j = 1:Nst
    W_line_h = [0 8]
    FoS_line_h = [FoS(j,idW) FoS(j,idW)]
    plot(W_line_h,FoS_line_h,"r--","linewidth",1)
    plot(8,FoS(j,idW),"ro","linewidth",4.5)
end

// allowable factor-of-safety for ductile material
fos_min_x = [0 12]; fos_min_y = [3 3]
plot(fos_min_x,fos_min_y,"k--","linewidth",1.2)

// annotation or labelling
eqn_1 ="$\sigma_d=\frac{F}{A},\;FoS=\frac{\sigma_u}{\sigma_d},\;FoS=\frac{\sigma_y}{\sigma_d}$"
ann_1 = xstring(6.25,19.0,eqn_1)
ann_1.font_size = 3.5
eqn_2 = "$\text{Based on} \; \sigma_u \, , \; FoS_{allow} \geq 3$"
ann_2 = xstring(6.75,17.0,eqn_2)
ann_2.font_size = 3.5
eqn_3 = "$\text{Based on} \; \sigma_y \, , \; FoS_{allow} \geq 1.5$"
ann_3 = xstring(6.75,15.0,eqn_3)
ann_3.font_size = 3.5

ann_4 = xstring(0,FoS(1,idW),"$FoS_{#1} = 4.1$")
ann_4.font_size = 3.5
ann_5 = xstring(0,FoS(2,idW),"$FoS_{#2} = 8.2$")
ann_5.font_size = 3.5
ann_6 = xstring(0,FoS(3,idW),"$FoS_{#3} = 12.3$")
ann_6.font_size = 3.5

ann_7 = xstring(0,fos_min_y(1)-1.75,"$FoS_{allow} = 3.0$")
ann_7.font_size = 3.5

Feel free to comment for inquiry, clarification, correction or suggestion for improvement. Drop your email to make a request to the author.

Disclaimer: The formulas and calculations presented are for technical reference only. Users must verify the accuracy and ensure compliance with applicable engineering standards, codes, and safety requirements before practical application.

References

1. Virgil Moring Faires. Design of Machine Elements. 4^th Ed. The Macmillan Company, New York. 1968.
2. Robert H. Cramer. Machine Design. 3^rd Ed. Addison-Wesley Publishing Company, Inc. 1984. 
3. F.L. Singer and A. Pytel. Strength of Materials. 3^rd Ed. Harper & Row, Publishers, New York. 1980. 
4. F.L. Singer. Engineering Mechanics. 2^nd Ed. Harper International Edition, New York, Evanston & London. 1970. 

A Shooting Method Approach to Boundary-Value ODEs

Topic: Differential Equations Boundary-Value Problems

Subject: Numerical Methods

Tool: Scilab

By: Gani Comia, July 2025

• 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 2^nd-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.

Disclaimer: The formulas and calculations presented are for technical reference only. Users must verify the accuracy and ensure compliance with applicable engineering standards, codes, and safety requirements before practical application.

References

1. Gerald, C. and Wheatly, P., Applied Numerical Analysis. Addison-Wesley Publishing Company, Inc. 1994.
2. ChatGPT. AI Tools. Accessed on July 12, 2025.
3. Google AI Overview. AI Tools. Accessed on July 12, 2025.