C++ and Python for Numerical Solution to Spring-Mass-Damper Model

Topic: SDOF Spring-Mass-Damper System
Subject: Numerical Methods & Mechanical Vibration
Tools: C++, CodeBlocks, GitHub Copilot, LabPlot, Python, VS Code, and QCAD

By: Gani Comia, May 2025

• Spring-Mass-Damper Model and Tools for Numerical Solution

Numerical solutions, such as Euler's method or Runge-Kutta method, can be used to determine the free response of a damped single-degree-of-freedom system. These frequently employed numerical techniques are utilized to solve vibration problems for which closed-form analytical solutions are not readily obtainable. Runge-Kutta simulations provide a higher degree of accuracy and quality compared to those generated by the Euler method.

C++ and Python programming language will be used in this article to solve numerically the response from a Spring-Mass-Damper Model using the Runge-Kutta method of 4^th-Order. One may refer to an article in this blog entitled “Simple Harmonic Motion Solution using Runge-Kutta RK4” for the discussion of RK4 numerical method as applied to 2^nd-Order ODE.

As observed, free oscillations in most systems eventually decay and the motion ceases entirely. This phenomenon can be modelled by adding a damper to the spring-mass model of the mechanical systems. The equation of motion of the spring-mass-damper model is represented by

$$m \, y''(t) + c \, y'(t) + k \, y(t) = 0  \tag{1}$$

Where:

\(m\) – mass, \(kg\)

\(c\) – damping coefficient, \(kg/s\)

\(k\) – spring constant, \(N/m\)

\(y’’(t)\) – acceleration as a function of time, \(m/s^2\)

\(y’(t)\) – velocity as a function of time, \(m/s\)

\(y(t)\) – displacement as a function of time, \(m\)

Subject to initial conditions

$$y(0) = y_0 \quad \text{and} \quad y'(0) = v_0  \tag{2}$$

Where:

\(y_0\) – initial position, \(m\), at \(t = 0 \; s\)

\(v_0\) – initial velocity, \(m/s\), at \(t = 0 \, s\)

Equation (1) is also referred to as a damped single-degree-of-freedom (SDOF) system. In analyzing the damped SDOF system, the non-dimensional number, \(\zeta\), called the damping ratio dictates the type of its response. Damping ratio is defined in Equation (3)

$$\zeta = \frac{c}{2 \, \sqrt{k \, m}}  \tag{3}$$

Damping ratio defines the types of motion of the response as follows

1. Underdamped: \(0 < \zeta < 1\)
2. Overdamped: \(\zeta > 1\)
3. Critically Damped:  \(\zeta = 1\)

• Application

Figure 1 shows a sample application where the response in terms of displacement and velocity as a function of time are to be solved.

Figure 1. Schematic of a SDOF Spring-Mass-Damper Vibration Model.

The given sample of damped system can be represented by the Equation (4).

$$3 \, y''(t) + y'(t) + 2 \, y(t) = 0  \tag{4}$$

With initial conditions as follows

$$y(0) = 0 \quad \text{and} \quad y'(0) = 0.25  \tag{5}$$

• Numerical Solution Using C++, Codeblocks, Github Copilot, and LabPlot.

Figure 2 is the plot of the numerical solution or the response on the system in terms of displacement and velocity. The code is a modified solution through the assistance of Github Copilot (AI). Codeblocks is used as the Integrated Development Environment (IDE) for the C++. The calculation results from C++ are exported as .csv file. LabPlot is the application used to import the .csv file to make the plot as shown in Figure 2.

Figure 2. Plot of Numerical Solution to Damped System using C++ and LabPlot.

• C++ Code

#include <iostream>
#include <fstream>
#include <vector>
#include <cmath>
#include <iomanip>

// Define the system of ODEs
// y1' = y2
// y2' = f(t, y1, y2) = -(2/3)*y1-(1/3)*y2

double f1(double t, double y1, double y2) {
    return y2; // y1' = y2
}
double f2(double t, double y1, double y2) {
    // Example: y'' = -(2/3)*y1-(1/3)*y2 (Simple Harmonic Oscillator)
    return -(2.0/3.0)*y1-(1.0/3.0)*y2; // y2' = -(2/3)*y1-(1/3)*y2
}
// Runge-Kutta 4th Order Method
void rungeKutta4(double t0, double y10, double y20, double tEnd, double h, const std::string& filename) {
    double t = t0;
    double y1 = y10;
    double y2 = y20;
    // Open the CSV file for writing
    std::ofstream file(filename);
    if (!file.is_open()) {
        std::cerr << "Error: Could not open file " << filename << " for writing.\n";
        return;
    }
    // Write the header to the CSV file
    file << "t,y1,y2\n";
    // Write the initial conditions
    file << t << "," << y1 << "," << y2 << "\n";
    while (t < tEnd) {
        // Compute RK4 coefficients for y1
        double k1_y1 = h * f1(t, y1, y2);
        double k1_y2 = h * f2(t, y1, y2);
        double k2_y1 = h * f1(t + h / 2.0, y1 + k1_y1 / 2.0, y2 + k1_y2 / 2.0);
        double k2_y2 = h * f2(t + h / 2.0, y1 + k1_y1 / 2.0, y2 + k1_y2 / 2.0);
        double k3_y1 = h * f1(t + h / 2.0, y1 + k2_y1 / 2.0, y2 + k2_y2 / 2.0);
        double k3_y2 = h * f2(t + h / 2.0, y1 + k2_y1 / 2.0, y2 + k2_y2 / 2.0);
        double k4_y1 = h * f1(t + h, y1 + k3_y1, y2 + k3_y2);
        double k4_y2 = h * f2(t + h, y1 + k3_y1, y2 + k3_y2);
        // Update y1 and y2
        y1 += (k1_y1 + 2 * k2_y1 + 2 * k3_y1 + k4_y1) / 6.0;
        y2 += (k1_y2 + 2 * k2_y2 + 2 * k3_y2 + k4_y2) / 6.0;
        // Update time
        t += h;
        // Write the results to the CSV file
        file << t << "," << y1 << "," << y2 << "\n";
    }
    // Close the file
    file.close();
    std::cout << "Results saved to " << filename << "\n";
}

int main() {
    // Initial conditions
    double t0 = 0.0;    // Initial time
    double y10 = 0.0;   // Initial value of y (y1)
    double y20 = 0.25;   // Initial value of y' (y2)
    double tEnd = 20.0; // End time
    double h = 0.1;     // Step size
    // Output file name
    std::string filename = "rk4_2nd_ode_output.csv";
    // Solve the ODE and save results to the CSV file
    rungeKutta4(t0, y10, y20, tEnd, h, filename);
    return 0;
}

• Numerical Solution Using Python, Github Copilot, and VS Code.

Figure 3 is the plot for the response of the same damped system also in terms of displacement and velocity. Github Copilot (Artificial Intelligence) is used to generate the Python script with minor modification. VSCode is used as IDE to run the script and plot the results.

Figure 3. Plot of Numerical Solution to Damped System using Python Script.

• Python Script

import numpy as np
import matplotlib.pyplot as plt

def runge_kutta_4_system(funcs, y0, t0, t_end, h):
    """

    Solve a system of first-order ODEs using the 4th-order Runge-Kutta method (RK4).
    Parameters:
        funcs (list): List of functions representing the system of ODEs.
        y0 (list): Initial values for the system [y1(0), y2(0), ...].
        t0 (float): Initial time.
        t_end (float): End time.
        h (float): Step size.
    Returns:
        tuple: Arrays of t and y values for each variable.
    """
    t_values = np.arange(t0, t_end + h, h)
    y_values = np.zeros((len(t_values), len(y0)))
    y_values[0] = y0
    for i in range(1, len(t_values)):
        t = t_values[i - 1]
        y = y_values[i - 1]
        k1 = [h * f(t, *y) for f in funcs]
        k2 = [h * f(t + h / 2, *(y + np.array(k1) / 2)) for f in funcs]
        k3 = [h * f(t + h / 2, *(y + np.array(k2) / 2)) for f in funcs]
        k4 = [h * f(t + h, *(y + np.array(k3))) for f in funcs]
        y_values[i] = y + (np.array(k1) + 2 * np.array(k2) + 2 * np.array(k3) + np.array(k4)) / 6
    return t_values, y_values

# Example application
if __name__ == "__main__":
    # Define the system of ODEs for 3y'' + y' + 2y = 0
    # Let y1 = y and y2 = -(2.0/3.0)*y1-(1.0/3.0)*y2, then:

    # y1' = y2
    # y2' = -(2.0/3.0)*y1-(1.0/3.0)*y2
    def f1(t, y1, y2):
        return y2
    def f2(t, y1, y2):
        return -(2.0/3.0)*y1-(1.0/3.0)*y2

    # Initial conditions
    y0 = [0, 0.25]  # y(0) = 0, y'(0) = 0.25
    t0 = 0       # Initial time
    t_end = 20   # End time
    h = 0.1      # Step size
    # Solve the system
    t_values, y_values = runge_kutta_4_system([f1, f2], y0, t0, t_end, h)
    # Extract y1 (the solution to the original ODE)
    y1_values = y_values[:, 0]
    y2_values = y_values[:, 1]

    # Plot the solution
    plt.figure(figsize=(8, 6))
    plt.plot(t_values, y1_values, linewidth=3, label="y(t)", color="blue")
    plt.plot(t_values, y2_values, linestyle='--', linewidth=2, label="y'(t)", color="red")
    plt.title("Using Python to Solve 2nd-Order ODE with RK4")
    plt.xlabel("Time, t")
    plt.ylabel("Displacement, y(t) & Velocity, y'(t)")
    plt.grid(True)
    plt.ylim(-0.3, 0.3)
    plt.legend()
    plt.show()

The damping ratio for the given application can be solved as shown

$$\zeta = \frac{c}{2 \, \sqrt{k \, m}} = \frac{1}{2 \sqrt{(2)(3)}} = 0.204$$

With the damping ratio, \(\zeta = 0.204\), the motion of the response of the system is considered underdamped.

C++ and Python programming languages using Codeblocks and VS Code as IDEs, and LabPlot for plotting application are some powerful engineering tools to analyze a system of mechanical vibration. AI through the use of Github Copilot provided the applicable codes for the solution.

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

References

1. D.I.Inman. Engineering Vibration. 2^nd Ed. Prentice Hall International, Inc., New Jersey 07458. 2001.
2. GitHub Copilot. AI Tools. Accessed on May 8, 2025 within the VS Code environment. 

Simple Harmonic Motion Solution using Runge-Kutta RK4 Method

Topic: Oscillatory Motion and Numerical Solution of ODE

Subject: Physics and Numerical Methods

Tool: Scilab

By: Gani Comia, Jan. 2025

The 4^th -order Runge-Kutta methods (RK4) is also used for numerical technique of approximating solution to the initial value problem (IVP) of 2^nd- order ordinary differential equation (ODE). This article presents the solution to the basic equation of simple harmonic motion (SHM) using RK4 algorithm in Scilab.

SHM is a form of nonhomogeneous linear differential equation (LDE) of second-order as shown in Equation (1).

$$y''\;+\;p(t)\,y'\;+\;q(t)\,y\;=\;r(t) \tag{1}$$ 

In order to solve a given higher-order differential equation, it is required to be converted to a system of first-order differential equations as presented below.

Let us define the following dependent variables, \(y_1\) and \(y_2\):

$$y_1\;=\;y \tag{2}$$

$$y_2\;=\;y' \tag{3}$$

Then, the nonhomogeneous LDE can be re-written as a system of 1^st-order ODE in Equation (4) and (5):

$$y'_1\;=\;y_2 \tag{4}$$

$$y'_2\;=\;-p(t)\,y_2\;-\;q(t)\,y_1\;+\;r(t) \tag{5}$$

With that presented, let us define and illustrate a problem using the basic equation of SHM and convert its 2^nd-order ODE to a system of 1^st-order to use the RK4 method in finding its approximate solution.

Figure 1. Single-Degree of Freedom (SDOF) Spring-Mass Model Oscillation Problem

The SDOF spring-mass model oscillation can be described by the basic equation of SHM shown in Equation (6).

$$m\,\frac{d^2y}{dt^2}\,+\,k\,y\,=\,0, \quad y(t_0)=y_0\;\;\text{and}\;\;y'(t_0)=y'_0 \tag{6}$$

Where:

\(m\) – mass

\(k\) – spring constant

\(t\) – time

\(t_0\) - initial time

\(y_0\) – initial displacement

\(y’_0\) – initial velocity

The solution to the ODE of SHM can be approximated using the Runge-Kutta (RK4) Method. The RK4 will be using the same algorithm presented in the previous article “First-Order IVP ODE Solution using Runge-Kutta RK4 Method” with the link shown: https://gani-mech-toolbox.blogspot.com/2025/01/first-order-ivp-ode-solution-using.html. Figure 2 shows the solution to the problem illustrated in Figure 1.

Figure 2. Solution to Spring-Mass Model System with RK4 using Scilab Script

The Scilab script to solve the given problem is presented below. Each block of the code is commented comprehensively for better understanding of the program sequence by the reader.

Scilab Script:

clear;clc;
// single-degree of freedom spring-mass model oscillator system
// m*y'' + k*y = 0  or  y'' + (k/m)*y = 0
// nonhomogeneous linear differential equation (LDE) of second-order
// y'' + p(t)y' + q(t)y = g(t) 
// where: p(t)=0; q(t)=k/m; g(t)=0

// (1) 2nd-order LDE re-written in the system of 1st-order ODEs
function dydt=f(t, y)
    m = 1;                  // kg, mass of the object
    k = 20;                 // N/m, spring constant
    p = 0;                  // coefficient of y'
    q = k/m;                // coefficient of y
    g = 0;                  // forcing function
    dydt = [y(2); -p*y(2) - q*y(1) + g];
endfunction

As you may compare the \(rk4\) function in block (2) to the previous article, the final time, \(t_f\), is used as one of the input arguments instead of the number of intervals, \(n\). Both \(rk4\) functions can be used to come up with the same solution.

// (2) RK4 function or algorithm for 1st-order ODEs
function [t, y]=rk4(f, t0, y0, h, tf)
    n = round((tf - t0)./h);
    t(1) = t0;          // initial time, t0
    y(:, 1) = y0;       // y(t0)=y0 & y'(t0)=y0' in column vector

    for i = 1:n
        // approximate values of slopes of solution curves
        k1 = h * f(t(i), y(:, i));
        k2 = h * f(t(i) + h / 2, y(:, i) + k1 / 2);
        k3 = h * f(t(i) + h / 2, y(:, i) + k2 / 2);
        k4 = h * f(t(i) + h, y(:, i) + k3);
        // estimates of numerical solution
        y(:, i + 1) = y(:, i) + (k1 + 2 * k2 + 2 * k3 + k4) / 6;
        t(i+1) = t(i) + h;
    end
endfunction

// (3) initial conditions and parameters
t0 = 0;             // sec, initial time
tf = 5;             // sec, final time
h = 0.01;           // sec, time step size
y0 = [0.05; 0];     // initial disp & velocity: y(0)=0.05, y'(0)=0

// (4) solution to the ODE by calling RK4 defined function
[t, y] = rk4(f, t0, y0, h, tf);

// (5) solution plot
clf;
plot(t, y(1, :),"b-","linewidth",4);       // displacement plot
plot(t, y(2, :),"r--","linewidth",1.8);    // velocity plot
xlabel("$\Large t \quad \text{(time)}$");
ylabel("$\Large y(t)\;\;\text{or}\;\;\frac{dy(t)}{dt}$");
note1 = "Spring-Mass Model System using Runge-Kutta RK4";
title(note1,"fontsize",3.5);
note2 = ["$\Large y(t)$","$\Large \frac{dy(t)}{dy}$"];
legend(note2, with_box = %F);
note3 = "https://gani-mech-toolbox.blogspot.com"
xstring(0.1,0.35,note3);
xgrid(3,1);

// (6) other plot information
disp(max(y(1,:)), min(y(1,:)))
disp(max(y(2,:)), min(y(2,:)))
ax = gca();
ax.data_bounds = [0 -0.4; 5 0.4];

// (7) natural frequency
format(5);
k = 20;                     // N/m, spring constant
m = 1;                      // kg, mass of object
omega = sqrt(k/m);          // Hz, natural frequency
disp(omega);
note4 = ["$\LARGE\mathbf{\omega_n}\;=\;\text{4.5 Hz}$"];
xstring(0.1,0.25,note4)

The Scilab script in this article can be one of your references to solve any engineering problems that can be described by the nonhomogeneous linear differential equation (LDE) of second-order.

Feel free to comment for inquiry, clarification, or suggestion for improvement. Drop your email to request the soft copy of the file.

First-Order IVP ODE Solution using Runge-Kutta RK4 Method

Topic: Numerical Solution to Ordinary Differential Equation

Subject: Numerical Methods

Tool: Scilab

By: Gani Comia, Jan. 2025

The 4^th -order Runge-Kutta methods (RK4) is a popular and most widely used numerical technique of approximating solution to the initial value problem (IVP) of 1^st -order ordinary differential equation (ODE). This article presents the equations and algorithm to apply the RK4 to a problem using Scilab scripts.

The RK4 method estimates the solution \(y=f(x)\) of an ODE with initial value condition (IVC) of the form:

$$\frac{dy}{dx}\,=\,f(x,y), \quad y(x_0)\,=\,y_0  \tag{1}$$

The RK4 algorithm requires a set of approximate values \((k_1, \; k_2, \; k_3, \; \text{and} \; k_4)\) of the slope of the solution curve at different points within the current step over an interval, \(h\).

$$\begin{aligned}k_1 &= h \cdot f \left( x_n, y_n \right) \\ k_2 &= h \cdot f \left( x_n + \frac{h}{2}, y_n+\frac{k_1}{2} \right)  \\ k_3 &= h \cdot f \left( x_n + \frac{h}{2}, y_n + \frac{k_2}{2} \right) \\ k_4 &= h \cdot f \left( x_n + h , y_n + k_3 \right) \end{aligned} \tag{2}$$

The \(k\) values from Equation (2) are then weighted and combined using Equation (3) to determine the estimates for \(y_{n + 1}\) of the numerical solution and \(x_{n + 1}\) is updated as per Equation (4).

$$y_{n+1} = y_n + \frac{1}{6} \left( k_1 + 2 k_2 + 2 k_3 + k_4 \right)  \tag{3}$$

$$x_{n+1} = x_n + h  \tag{4}$$

Let us apply using Scilab scripting the implementation of approximating solution to the ODE with initial condition in Equation (5) using the RK4 solution algorithm.

$$\frac{dy}{dx} = \frac{2x}{y} - xy , \quad y(0) = 1  \tag{5}$$

• Scilab Script

The first block of the script is the function definition in Scilab of Equation (5).

clear;clc;

// (1) 1st-order ODE (dy/dx = f(x, y)) in Scilab form
function dydx=f(x, y)
    dydx = (2*x./y)-(x.*y);     // Example: dy/dx = 2x/y - xy
endfunction

The RK4 solution algorithm is defined in Scilab and shown in the second block of script.

// (2) Scilab function of Runge-Kutta 4th order method (RK4)
function [x, y]=rk4(f, x0, y0, h, n)
    x = zeros(1, n+1);
    y = zeros(1, n+1);
    x(1) = x0;
    y(1) = y0;

    for i = 1:n
        k1 = h * f(x(i), y(i));
        k2 = h * f(x(i) + h/2, y(i) + k1/2);
        k3 = h * f(x(i) + h/2, y(i) + k2/2);
        k4 = h * f(x(i) + h, y(i) + k3);

        x(i+1) = x(i) + h;
        y(i+1) = y(i) + (k1 + 2*k2 + 2*k3 + k4) / 6;
    end
endfunction

Initial conditions and RK4 parameters such as step size and number of intervals are shown below.

// (3) initial conditions and parameters
x0 = 0;                         // initial value of x
y0 = 1;                         // initial value of y
h = 0.1;                        // step size
n = 50;                         // no. of interval

Scilab function for the RK4 algorithm is executed to find the solution vectors for both \(x\) and \(y\) variables. Displaying on the Scilab console the last element of the variables provides an insight for the domain and range of the graph.

// (4) numerical solution to the 1st-order ODE IVP using RK4
[x, y] = rk4(f, x0, y0, h, n); 
disp(x($),y($));        // console display of the last value of x & y

Finally, visualization of numerical solutions requires plotting values of variables using appropriate graph. Shown in Figure 1 is the solution to the given ODE with IVC at \((0,1)\).

// (5) plot of the solution
clf;
f = gcf();
f.figure_size = [600,600];

plot(x, y,"b-","linewidth",3);
xlabel("$\large\textbf{x-value}$");
ylabel("$\large\textbf{y-value}$");
title("$\large\text{Solution to}\;\frac{dy}{dx}=\frac{2x}{y}-{xy}$");
legend("$\LARGE{y\,=\,f(x)}\;\text{with RK4}$",with_box=%F);
xstring(2.75,1.05,"https://gani-mech-toolbox.blogspot.com");
xgrid(3,1);

ax = gca();
ax.data_bounds = [0 1; 5 1.5];

Scilab Output (Figure 1)

Figure 1. Solution Plot of 1st-order ODE with IVC using RK4.

RK4 solution to 2^nd-order ODE with IVC shall be presented in our next article.

Feel free to comment for inquiry, clarification, or suggestion for improvement. Drop your email to request the soft copy of the file.

Newton’s Law of Cooling with Euler’s Method Solution to IVP ODE

Topic: Newton’s Law of Cooling and Euler’s Method to IVP ODE

Subject: Physics and Numerical Methods

Tool: Scilab

By: Gani Comia, Jan. 2025

This article presents a solution as to how long it will take for an object to cool down to a specified temperature and thereby to ambient temperature. The object is being taken out of the oven at \(230^\circ C\) and put in place where the ambient temperature is \(30^\circ C\). Neglecting the heat transfer due to conduction where the object is in contact with, the time to cool it down to a particular temperature due to natural convection and radiation can be estimated using the ordinary differential equation (ODE) representing the Newton’s Law of Cooling shown in Equation (1).

$$\frac{dT}{dt}\;=\;-k \left(T(t) - T_{env} \right) \tag{1}$$

Where:

\(T(t)\) - object temperature at time \(t\),

\(T_{env}\) - environment temperature,

\(k\) - cooling constant, 

\(dT/dt\) - rate of change of temperature with time,

The analytical solution to the ODE representing Newton’s law of cooling can be found by integrating both sides of the equation as shown in Equation (2).

$$\int_{T_{0}}^{T_{t}} dT\;=\;-k \int_{0}^{t} \left(T(t)-T_{env} \right) dt \tag{2}$$

For this article, a numerical solution using the Euler method of approximating the solution of a first-order differential equation is presented with the use of Scilab. Euler methods require an ODE to be in the form of Equation (3).

$$\frac{dT}{dt}\;=\;f\left(t,T\right) \tag{3}$$

With the given initial condition,

$$\left(t_{0},\;T_{0}\right) \tag{4}$$

The approximate solution using Euler's method is shown in Equation (5).

$$\frac{\Delta T}{\Delta t}\; \approx \; \frac{T_{n+1}-T_{n}}{h}\; \approx \; f\left(t,T \right) \tag{5}$$

Or by rearranging and solving for \(T_{n+1}\),

$$T_{n+1}\;=\;T_{n}+h \;f\left(t,T\right) \tag{6}$$

The \( f(t,T)\) requires finding the cooling constant, \(k\), for a specific material. If this property seems to be unavailable or difficult to find, experimental data can be used as the basis for calculating the \(k\) using Equation (7).

$$k\;=\;- \frac{1}{t} \ln \left[ \frac{T(t) - T_{env}}{T_{0} - T_{env}} \right] \tag{7}$$

Figure (1) is the illustration of the condition to be given a solution using Newton’s law of cooling by applying Euler's method of solving ODE.

Figure 1. An object at \(230^\circ C\) exposed in an environment at \(30^\circ C\).

To provide a simulation of the given condition in Figure (1), experimental data were gathered to calculate the cooling constant, \(k\), from the given temperature and time. The cooling constant, \(k = 0.08 \;\text{per min}\), was calculated from the given parameters and Scilab script.

// Copyright (C) 2024 - Gani Comia
// Experimental Cooling Constant
clear;clc;

// (1) defining cooling constant function
function k=f(t, Tt, T0, Tenv)
    // where:
    // t - min, time consideration
    // Tt - deg C, temperature at time t
    // T0 - deg C, initial temperature
    // Tenv - deg C, environment temperature
    k = -(1/t).*log((Tt - Tenv)./(T0 - Tenv));
endfunction

// (2) experimental data
t = 30;                     // min, time between T0 and Tt
T0 = 230;                   // deg C, initial temp
Tt = 50;                    // deg C, temp at time t
Tenv = 30;                  // deg C, ambient temp

// (3) cooling constant calculation
k = f(t,Tt,T0,Tenv)
mprintf("Cooling Constant, k = %4.2f per min\n", k)

Having this cooling constant, a simulation of Newton’s law of cooling can be performed. The codes provided in the simulation should be saved on a different file from the script used for estimating the cooling constant.

The first block of the script defines the physical parameter under consideration.

clear;clc;

// (1) physical parameters
T0 = 230;                   // deg C, initial temp of the object
Tenv = 30;                  // deg C, ambient temp
k = 0.08;                   // per min, experimental cooling constant
dt = 0.05;                  // min, time step
totalTime = 80;             // min, total simulation time

The time domain in the simulation is divided into smaller step sizes for accuracy of solution. Initial time and temperature are given on the next block of script.

// (2) initial conditions and derived parameters
time = 0:dt:totalTime;
nSteps = length(time);
T = zeros(nSteps, 1);
T(1) = T0;

Equation (6), which is the numerical solution to the ODE using the Euler method, is transformed in the Scilab script. It uses the \(for-loop\) function to obtain the solution over the desired range.

// (3) numerical solution to ODE using Euler Method
for i = 1:nSteps-1
    dTdt = -k * (T(i) - Tenv);  // deg C/min, rate of temp change
    T(i+1) = T(i) + dTdt * dt;  // deg C, temp at the next time step
end

The Scilab \(plot()\) function is used to visualize the solution from the Euler method of approximation. Finding the temperature for a given time can be estimated using interpolation by the use of the \(interp1()\) function.

// (4) graph for visualization of solution
clf;
plot(time, T, "b","linewidth",2.5);
xlabel("Time, t [ min ]","fontsize",2);
ylabel("Temperature, T [ °C ]","fontsize",2);
title("Newton''s Law of Cooling Simulation","fontsize",3);
note = ["T0="+string(T0)+"°C, Tenv="+string(Tenv)+"°C, k="+string(k)]
legend(note,with_box = %F);
xgrid(3,1);

// (5) interpolating temp, Tx, for the particular time, tx
format(6)
tx = 17.5;                      // min, sample time
Tx = interp1(time, T, tx);      // deg C, interpolated temp
disp(tx, Tx)

// (6) plotting lines for tx and Tx
xptVer = [tx tx];
yptVer = [0 Tx];
plot(xptVer, yptVer, "r--","linewidth",1.5)
xptHor = [0 tx];
yptHor = [Tx Tx];
plot(xptHor, yptHor, "r--","linewidth",1.5)
plot(tx,Tx,"ro");
note1 = ["t (min) ="+string(tx)+"; T (°C) ="+string(Tx)];
xstring(tx,Tx,note1);
xstring(45,200,["https://gani-mech-toolbox.blogspot.com"])
eqn1=["$\LARGE\frac{dT}{dt}\;=\;{-k}\left(T-T_\mathbf{{env}}\right)$"];
xstring(45,166,eqn1);
a = gca();
a.children.children.mark_background = 9; 

Figure 2 is the solution over the range of time.

Figure 2. Plot of the Solution to Newton’s Law of Cooling using Euler’s Method.

This Scilab toolbox is a powerful tool for visualizing solution data and calculating specific conditions.

Feel free to comment for inquiry, clarification, or suggestion for improvement. Drop your email to request the soft copy of the file.

Laplace Equation 1D for Steady-State Heat Conduction Model in BVP PDE

Topic: Steady-State Heat Conduction

Subject: Heat Transfer and Numerical Methods

Tool: Scilab

By: Gani Comia, Jan. 2025

The steady-state heat conduction model, as represented by the Laplace equation in one dimension (1D), is the simplest partial differential equation (PDE) you may encounter solving. This equation explains the temperature characteristics or profile between the boundaries of the given domain. Presented in this article is the application of the Laplace equation in 1D using the Finite Difference Method (FDM) for the solution.

The Laplace equation in 1D for the steady-state heat conduction is shown in Equation (1).

$$\frac{\partial^2 u}{\partial x^2}\;=\;0 \tag{1}$$

The second-order PDE in 1D can also be written in ordinary differential equations (ODE). This can be approximated using the FDM shown in Equation (2).

$$\frac{d^2 u}{d x^2}\; \approx \; \frac{u_{i+1} - 2 u_i + u_{i-1}}{\left(\Delta x\right)^2}\;=\;0 \tag{2}$$

This leads to the solution in Equation (3),

$$u_i\;=\; \frac{u_{i+1} + u_{i-1}}{2} \quad for \quad i\;=\;2,3, \ldots ,N-1 \tag{3}$$

With boundary conditions,

$$u_1\;=\;u_1, \quad u_N\;=\;u_L \tag{4}$$

Let us illustrate an example of steady-state heat conduction as shown in Figure (1) and apply the FDM solution to understand its temperature profile between boundaries. The Scilab script to implement the solution will be explained in detail in relation to the Laplace equation in 1D.

Figure 1. Sample Domain under Steady-State Heat Conduction

The solution using the Scilab script is shown below. A short explanation is provided for each code snippet.

Scilab Script.

clear;clc;clf;

// (1) primary & derived parameters
L = 100;                            // mm, length of workpiece
Nx = 11;                            // unit, number of nodes
dx = L / (Nx - 1);                  // mm, spatial step

Equation (5) is the boundary condition for this problem.

$$u_1\;=\;u_1, \quad u_{2,3, \dots ,N-1,N}\;=\;u_L \tag{5}$$

// (2) boundary temperature conditions
u = zeros(Nx, 1) + 25;              // node temp cond @ 25 deg C
u(1) = 100;                         // boundary temperature

The temperature condition at node points between boundaries can be solved using the "for-end” function of the Scilab.

// (3) temperature profile at node points using FDM
for i = 2:Nx-1
    u(i) = (u(i+1) + u(i-1))/2;
end

Below is the script for plotting the temperature profile.

// (4) plotting of temperature at node points
f = gcf()
f.figure_size = [600,600]

len = linspace(0,L,Nx);
plot(len, u, 'r-o');
title("Temperature Profile @ Eleven (11) Node Points","fontsize",3)
xlabel("Length, x [ mm ]","fontsize",3)
ylabel("Temperature, T [ C ]","fontsize",3)
xstring(58,73,"https://gani-mech-toolbox.blogspot.com")
eqn = ["$\LARGE \;\frac{\partial^2\,T}{\partial x^2}\,=\,0$"]
legend(eqn,with_box=%f)
xgrid(3,1)

a = gca()
a.data_bounds = [0 0; L max(u)]
a.children.children.mark_background = 9;

Figure (2) is the resulting plot from executing the whole script.

Figure 2. Graphical Solution for Laplace 1D using FDM with 11 Node Points

The given boundary conditions in Equation (5) will result in a curved-like temperature profile. If Equation (4) is considered as boundary conditions with the corresponding change in primary and derived parameters, specifically reducing the number of node points to 3, a linear temperature profile will become the solution.

Below is the modified script for 3 node points.

// (1) primary & derived parameters
L = 100;                            // mm, length of workpiece
Nx = 3;                             // unit, number of nodes
dx = L / (Nx - 1);                  // mm, spatial step

Figures (3) and (4) are the illustration and graphical solution for 3 node points leading to a linear temperature profile.

Figure 3. Sample Domain with 3 Node Points under Steady-State Heat Conduction

Figure 4. Graphical Solution to Laplace 1D using FDM for 3 Node Points

The number of node points taken into consideration in solving the Laplace equation using FDM dictates the type of temperature profile, linear or curve-like. You may send a comment about the consideration for a basis as to how many node points are to be used for certain conditions of the steady-state heat conduction model.

Feel free to comment for inquiry, clarification, or suggestion for improvement. Drop your email to request the soft copy of the file.