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.

Exponential Illusions: Modeling Population Growth with Caution

Topic: Exponential Growth Model and Limitation

Subject: Numerical Methods

Tool: Scilab

By: Gani Comia, July 2025

• Ordinary Differential Equation (ODE) for Exponential Growth

The ODE representing exponential growth is a fundamental model frequently used in real-world application of predictive modeling, particularly in demographics, ecology, and urban planning, to forecast future population sizes based on past trends and other factors.

The population growth is represented by the differential equation shown in Equation (1). The equation signifies that the rate of population change is directly proportional to the current population size.

$$\frac{dP}{dt}=r \, P \tag{1}$$

Where:

\(dP/dt\) – rate of change of the population with respect to time

\(P\) – population size at any given initial time

\(r\) – growth rate, an assumed positive constant that influences population growth

Equation (1) has a general solution as shown in Equation (2).

$$P(t) = P_0 \; e^{r \, t}, \quad P(0) = P_0 \tag{2}$$

Where:

\(P(t)\) – population at any given time, \(t\)

\(P_0\) – population at time, \(t = 0\)

The ODE representing the population growth model has the following limitations if applied in predictive modelling. First, it predicts unlimited and unrealistic growth over the long period of time, which is not sustainable in the long run due to factors resulting to limitation of resources. Second, the growth rate considered constant does not account for changes in birth and death rates or other environmental factors.

Below is the list of models or techniques for estimating or describing the population aside from exponential growth.

1.  Exponential Growth
2. Logistic Growth
3. System Dynamics
4. Naive Bayes and Decision Tree
5. Grey Prediction

• Numerical Solution to ODE using Euler's Method

The Scilab scripts used for calculations implements the Euler method to numerically approximate the solution of a first-order ordinary differential equation modeling exponential growth. Euler method requires an ODE and initial condition in the form of Equation (3).

$$\frac{dP}{dt} = f(t,P) , \quad (t_0 , P_0) \tag{3}$$

The approximate solution using Euler’s method is shown in Equation (4).

$$\frac{\Delta P}{\Delta t} \approx \frac{P_{n+1} - P_n}{h} \approx f(t, P) \tag{4}$$

Where:

\(h\) – step size, \(h = \Delta t\)

\(P_{n+1}\) is the solution to the ODE and can be solved by rearranging the Equation (4) leading to Equation (5).

$$P_{n+1} = P_n + h \; f(t, P) \tag{5}$$

The \(f(t,P)\) requires finding the growth rate, \(r\), for the two conditions as shown in Equation (6).

$$\large r = \frac{\ln \left( {\frac{P_t}{P_0}} \right)}{t} \tag{6}$$

• Philippine Statistics Authority (PSA) Population Estimates

This article will use the exponential growth model to predict and estimates the population growth of the Philippines in comparison with a country’s agency (PSA) estimates. The agency is mandated to conduct national censuses and surveys including those on population.

Table 1 presents the population projections from the PSA up to the year 2055. These projections will be analyzed alongside the model based on exponential growth.

Table 1. Philippine Population Projection until 2055

• PSA Estimates and Exponential Growth Comparison

Figure 1 illustrates the difference between two methods of estimation. The exponential growth model deviates from the PSA estimates after 25 years. Exponential growth models can accurately predict population trends during the initial growth phase, but their reliability decreases as other factors cause deviations from exponential behavior.

Figure 1. Population Growth Model and PSA Estimates

Below is the Scilab script to produce the results shown in Figure 1.

• Scilab Script

// Copyright (C) 2025 - I.S.Comia
// Date of creation: Jul 6, 2025
// Predictive Modeling using Exponential Growth
clear;clc;

// Subroutine
// Population or Exponential Growth Model (1st-Order ODE)
function dPdt=f(t, P)
    global r;
    dPdt = r * P;
endfunction

function r=growthRate(P_0, P_t, t)
    // calculate the estimated growth rate
    // input:
    //      P_0 - quantity at t = 0
    //      P_t - quantity at a given time t
    //      t - time
    // output:
    //      r - estimated growth rate
    r = log(P_t./P_0)./t;
endfunction

// Main function
// Problem case scenario
// Phil. Statistic Authority (PSA) estimated population from year 2020 to 2050
// Primary parameters
actYear = [2020 2025 2030 2035 2040 2045 2050 2055];
// Population in millions
Pop = [109.2027 113.8631 118.8738 123.9636 128.8260 133.0245 136.2989 138.6727];
Yrs = actYear - 2020;                       // years covered
t = Yrs(2)-Yrs(1);
global r;
r = growthRate(Pop(1),Pop(2),t);            // growth rate estimate
T_initial = Yrs(1);                         // initial time
P_initial = Pop(1);                         // initial population (millions)

// Secondary parameters
T_final = Yrs($);                           // years, final time
dt = 1/12;                                  // month, step size,
N = T_final/dt;                             // number of steps

// Euler method solution to IVP ODE
// Initial conditions for time and P-value
T(1) = T_initial;
P(1) = P_initial;
// Solution
for k = 1:N
    dPdt = f(T(k), P(k));
    T(k+1) = T(k) + dt;
    P(k+1) = P(k) + dPdt * dt;
end

// Plotting of ODE solution
clf;
fig = gcf()
fig.figure_size = [700, 700]
plot(T,P,"b-","LineWidth",4)
title("Population Growth Model, Year 2020~2055","FontSize",4)
xlabel("Time (years)","FontSize",3.5)
ylabel("Population (millions)","FontSize",3.5)
xgrid(7,1)

// Population estimate for T_value~5 years or Year 2025
idT = max(find(T <= 5.05));
T_value = T(idT);
P_value = P(idT);
mprintf("For T_value: %4.2f, P_value: %6.3f \n", T_value, P_value)

// Plot of information for Year 2025
plot(T_value,P_value,"rs","MarkerFaceColor","r","MarkerSize",13)
ann_1=xstring(T_value+0.75,P_value-1,["Year:2025,Pop:",string(round(P_value)),"M"])
ann_1.font_size = 2
ann_2 = xstring(20,102.5,"https://gani-mech-toolbox.blogspot.com")
ann_2.font_size = 2
ann_3 = xstring(T_initial-2.5,P_initial-2.5,"Year: 2020")
ann_3.font_size = 2

T_line = [T_value T_value T_initial-5];
P_line = [100 P_value P_value];
plot(T_line,P_line,"r--","LineWidth",1)

// PSA population estimate
plot(Yrs,Pop,"ko","MarkerFaceColor","k","MarkerSize",8)
ann_4 = ["Exponential Growth Model";"Est Pop @ Year 2025";"";"PSA Pop Estimate"]
ann_5 = legend(ann_4,2,with_box=%F)
ann_5.font_size = 2

ax = gca()
ax.data_bounds = [-5 100; 40 150]

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. James, G., et al., Modern Engineering Mathematics. Addison-Wesley Publishing Company, Inc. 1993.
2. “Philippine Population is Projected to be around 138.67 Million by 2055 under Scenario 2”. Philippine Statistics Authority. January 31, 2024. https://psa.gov.ph/content/philippine-population-projected-be-around-13867-million-2055-under-scenario-2
3. ChatGPT. AI Tools. Accessed on July 6, 2025.
4. Google AI Overview. AI Tools. Accessed on July 9, 2025.