Advanced Engineering Applications

In this chapter we present pratical engineering problems solved with combinations of concept discussed under linear algebra, optimization, ordinary differential equations, numerical laplace inversion and special functions

Dimensionless Water Influx Estimation

Water influx in an oil reservoir is the migration of water from an aquifer into the pore spaces of the reservoir rock containing oil. This water movement is primarily driven by pressure differences between the aquifer and the reservoir as the oil is produced and reservoir pressure declines. The water influx can provide pressure support, helping to maintain reservoir pressure and sustain oil production. Hence, understanding and accurate estimation of water influx is crucial for optimizing oil recovery strategies and the long-term economic viability of an oil field.

For use in material balance computation in edge drive configuration, reservoir engneering books provide plots for Wd as a function of dimensionless radius and time

For \(Rd \leq 4\)

For \(5 \leq Rd \leq 10\)

In an edge drive configuration with the aquifer closed at its outer boundary, the governing equation gives:

\[\cfrac{\partial P}{\partial t} = \cfrac{1}{r}\cfrac{\partial}{\partial r}\left(r \cfrac{\partial P}{\partial r} \right)\]

\[P(t = 0, r) = 0, P(t, r = 1) = 1, \cfrac{\partial P}{\partial r}(t, r = r_D) = 0\]

The solution in laplace space:

\[P(s, r) = \Phi_1 I_0(r\sqrt{s}) + \Phi_2 K_0(r\sqrt{s})\]

Using the boundary conditions to evaluate the constants and substitute them:

\[P(s, r) = \cfrac{K_1(r_D\sqrt{s}) I_0(r\sqrt{s}) + I_1(r_D\sqrt{s}) K_0(r\sqrt{s})}{s(K_1(r_D\sqrt{s}) I_0(\sqrt{s}) + I_1(r_D\sqrt{s}) K_0(\sqrt{s}))}\]

From Darcy law, we know that the rate of water influx is proportional to the negative rate of change of pressure with respect to radial position at the reservoir aquifer boundary, hence total water influx after a time t is thus:

\[W(t) = \int_{0}^{t_D} -\cfrac{\partial P}{\partial r} (\tau, r = 1) \partial \tau\]

This can be accomplised by performing the integration in laplace space before inverting to time space.

\[W(t) = \mathcal{L}^{-1}\left(\frac{-1}{s} \cfrac{\partial P}{\partial r} (s, r = 1) \right)\]

\[W(t) = \mathcal{L}^{-1}\left(\frac{1}{s\sqrt{s}} \cfrac{I_1(r_D\sqrt{s}) K_1(\sqrt{s}) - K_1(r_D\sqrt{s}) I_1(\sqrt{s})}{(I_1(r_D\sqrt{s}) K_0(\sqrt{s}) + K_1(r_D\sqrt{s}) I_0(\sqrt{s}))} \right)\]

Lets see how to compute water influx, and generate the started water influx plot as shown above

Note

The matlab implementation uses the numerical laplace inversion code from matlab file exchange[3]

CCL-MathPythonMatlab

CCL-Math Implementation

// import libraries
using System;
using System.Collections.Generic;
using SepalSolver;
using static SepalSolver.Math;

// define Wd function in time space.
double EdgeClosedBoundaryRadial_Wd(double tD, double rD)
{
    // define the embedded laplace space solution
    Func<double, double> lapW = new(s =>
    {
        double sqrts, sqrts3, rDsqrts, Num, Den;
        sqrts = Sqrt(s); sqrts3 = s * sqrts; rDsqrts = rD * sqrts;
        Num = BesselI(1, rDsqrts) * BesselK(1, sqrts) - BesselK(1, rDsqrts) * BesselI(1, sqrts);
        Den = BesselI(1, rDsqrts) * BesselK(0, sqrts) + BesselK(1, rDsqrts) * BesselI(0, sqrts);
        if (double.IsInfinity(Num) || double.IsInfinity(Den))
            return 1 / sqrts3;
        else
            return Num / (sqrts3 * Den);
    });
    return tD == 0 || rD == 1 ? 0 : NiLaplace(lapW, tD);
}

// define the time and radial mesh
double[] Rd = [2, 2.5, 3, 3.5, 4, 50];
ColVec Td = Logspace(-1, 2), Wd;
int end = Rd.Length - 1;
// compute the water influx and plot
hold = true;
List<string> lgd = [];
foreach (double rD in Rd)
{
    Wd = Td.Select(tD => EdgeClosedBoundaryRadial_Wd(tD, rD)).ToList();
    SemiLogx(Td, Wd, Linewidth: 2); lgd.Add("rD = " + rD);
}
lgd[end] = "rD = ∞";
Xlabel("tD"); Ylabel("WD");
Legend(lgd, Alignment.UpperLeft);
Axis([0.1, 100, 1, 8]);
Title("Dimensionless Water Influx");
SaveAs("Dimensionless-Water-Influx-CCL-Math-4dn.png");

// define the time and radial mesh
Rd = [5, 6, 7, 8, 9, 10, 50];
Td = Logspace(0, 3); end = Rd.Length - 1;

// compute the water influx and plot
hold = true;
lgd = [];
foreach (double rD in Rd)
{
    Wd = Td.Select(tD => EdgeClosedBoundaryRadial_Wd(tD, rD)).ToList();
    SemiLogx(Td, Wd, Linewidth: 2); lgd.Add("rD = " + rD);
}
lgd[end] = "rD = ∞";
Xlabel("tD"); Ylabel("WD");
Legend(lgd, Alignment.UpperLeft);
Axis([1, 1000, 0, 70]);
Title("Dimensionless Water Influx");
SaveAs("Dimensionless-Water-Influx-CCL-Math-5up.png");

Compressibility of Natural Gas

Natural gas compressibility is a measure of how much the volume of a given amount of natural gas decreases when pressure is applied to it. It is an important property to understand because it affects how much gas can be stored in a given space, and how easily the gas can be transported through pipelines.

Why Compressibility Matters

• Storage: Natural gas is often stored in underground reservoirs or tanks. The compressibility of the gas determines how much gas can be packed into these storage spaces at a given pressure.

• Transportation: In pipelines, higher pressure helps move gas efficiently. Compressibility affects how much pressure is needed to transport a certain amount of gas.

• Metering: Accurate measurement of natural gas flow relies on understanding its compressibility, as the volume of gas changes with pressure and temperature.

• Reservoir Engineering: In oil and gas reservoirs, compressibility is crucial for estimating gas reserves and predicting how the reservoir will perform under production.

Factors Affecting Compressibility

• Pressure: As pressure increases, the volume of natural gas decreases, meaning compressibility increases.

• Temperature: Higher temperatures generally decrease compressibility, as the gas molecules have more energy and resist being forced closer together.

• Gas Composition: Natural gas is a mixture of hydrocarbons, primarily methane. The exact composition affects compressibility, as different gases have different compressibility factors.

Calculating Compressibility

• Ideal Gas Law: At low pressures and high temperatures, natural gas behaves close to an ideal gas, and its compressibility can be estimated using the ideal gas law.

• Real Gas Equations: At higher pressures and lower temperatures, deviations from ideal behavior become significant. More complex equations, like the Peng-Robinson or Soave-Redlich-Kwong equations, are used to account for real gas behavior.

• Compressibility Factor (Z-factor): This factor corrects the ideal gas law for real gas behavior. It is often determined using charts or correlations based on pressure, temperature, and gas composition.

Importance of Accurate Compressibility Calculations

• Safety: Inaccurate compressibility calculations can lead to over-pressurization of storage vessels or pipelines, posing safety risks.

• Economic Efficiency: Optimizing pipeline and storage design, as well as accurately metering gas flow, requires precise compressibility data to ensure cost-effectiveness.

• Environmental Impact: Understanding gas behavior is essential for preventing leaks and ensuring the safe handling of natural gas, reducing its environmental impact.

Additional Points

• Isothermal Compressibility: This refers to compressibility measured at a constant temperature.

• Compressibility of Gas Mixtures: When dealing with natural gas, which is a mixture, the compressibility of the mixture needs to be determined based on the individual components and their proportions.

Z-factor

Here we show how to compute Z factor from Hall and Yarborough Ea derivative of Starline Carnahan’s Equation of state:

\[\begin{split}\begin{array}{rcl} t &=& 1/T_{pr} \\ A &=& 0.06125t\exp(-1.2(1-t)^2) \\ B &=& 14.76t - 9.76t^2 + 4.58t^3 \\ C &=& 90.7t - 242.2t^2 + 42.4t^3 \\ D &=& 2.18 + 2.82t \end{array}\end{split}\]

We solve for the root of:

\[-AP_{pr} + \frac{y+y^2+y^3-y^4}{(1-y)^3} - By^2+Cy^D = 0\]

And then compute:

\[z = \frac{AP_{pr}}{y}\]

CCL-MathPythonMatlab

CCL-Math Implementation

// import libraries
using System;
using System.Collections.Generic;
using SepalSolver;
using static SepalSolver.Math;

// define zfunction
static double ZfactorHY(double Pr, double Tr)
{
    // define variables
    double z = 1, t, tm1, tm1e2, t2, t3, A, B,
        C, D, r, y2, y3, y4, Den;

    // avoid computing z when Pr = 0.
    if (Pr != 0)
    {
        t = 1 / Tr; t2 = t * t; t3 = t2 * t;
        tm1 = 1 - t; tm1e2 = tm1 * tm1;
        A = 0.06125 * t * Exp(-1.2 * tm1e2);
        B = 14.76 * t - 9.76 * t2 + 4.58 * t3;
        C = 90.7 * t - 242.2 * t2 + 42.4 * t3;
        D = 2.18 + 2.82 * t; r = A * Pr;

        // define density equation
        var yfunc = new Func<double, double>(y =>
        {
            y2 = y * y; y3 = y2 * y; y4 = y3 * y;
            Den = Pow(1 - y, 3);
            return -A * Pr + (y + y2 + y3 - y4) / Den -
            B * y2 + C * Pow(y, D);
        });

        // solve the density equation
        var opts = SolverSet (StepFactor: 0.5 );
        double y = Fsolve(yfunc, r, opts);

        // compute the z factor
        z = A * Pr / y;
    }
    return z;
}

// set up ressure and temperature mesh
double[] Pr = Enumerable.Range(0, 301).Select(i => i * 0.05).ToArray();
double[] Tr = [1.05,    1.08,   1.12,   1.18,   1.26,   1.35,   1.47,
               1.61,    1.75,   1.91,   2.09,   2.29,   2.62,   3.00];
ColVec Z;

// compute z factors and plot them
List<string> Tlabels = [];
hold = true;
foreach( var tr in Tr)
{
    Z = Pr.Select(p => ZfactorHY(p, tr)).ToArray();
    Plot(Pr, Z);
    Tlabels.Add("Tr = " + tr);
}
Legend(Tlabels, Alignment.LowerRight );
SaveAs("Zfactor-Hall-Yarborough-CCL-Math.png");

Reduced Compressiblity

Reduced compressibility is given by:

\[c_rT_r = \frac{1}{Ppr} - \frac{1}{z}\frac{dz}{dPpr}\]

Using Hall and Yarborough Correlation, we can evaluate the reduced compressibility as:

\[\chi = \cfrac{1 + 4y + 4y^2 - 4y^3 + y^4}{(1 - y)^4} - 2By + CDy^{(D - 1)}\]

\[c_rT_r = \cfrac{A Tr}{y\chi};\]

CCL-MathPythonMatlab

CCL-Math Implementation

// import libraries
using System;
using System.Collections.Generic;
using SepalSolver;
using static SepalSolver.Math;

// define CrTr function
static double CrTrHY(double Pr, double Tr)
{
    double t, tm1, tm1e2, t2, t3, A, B, C, D, r, y2, y3, y4, yDm1, ym1p3, ym1p4, Den;

    // compute itnermediate variables
    t = 1 / Tr; t2 = t * t; t3 = t2 * t; tm1 = 1 - t; tm1e2 = tm1 * tm1;
    A = 0.06125 * t * Exp(-1.2 * tm1e2); B = 14.76 * t - 9.76 * t2 + 4.58 * t3;
    C = 90.7 * t - 242.2 * t2 + 42.4 * t3;  D = 2.18 + 2.82 * t; r = A * Pr;

    // define density function
    var yfunc = new Func<double, double>(y =>
    {
        y2 = y * y; y3 = y2 * y; y4 = y3 * y; ym1p3 = Pow(1 - y, 3);
        return -A * Pr + (y + y2 + y3 - y4) / ym1p3 - B * y2 + C * Pow(y, D);
    });

    // solve density function
    var opts = SolverSet(StepFactor: 0.5);
    if (Pr > 15) r *= 0.5;
    double y = Fsolve(yfunc, r, opts);

    // compute crtr
    y2 = y * y; y3 = y2 * y; y4 = y3 * y; ym1p4 = Pow(1 - y, 4); yDm1 = Pow(y, D - 1);
    Den = (1 + 4 * y + 4 * y2 - 4 * y3 + y4) / ym1p4 - 2 * B * y + C * D * yDm1;
    return A * Tr / (y * Den);
}

// set up ressure and temperature mesh
double[] Pr = Linspace(0.2, 25, 501);
double[] Tr = [1.05, 1.10, 1.20, 1.30, 1.40];
ColVec Cr;

// compute CrTr and plot them
List<string> Tlabels = [];
hold = true;
foreach( var tr in Tr)
{
    Cr = Pr.Select(p => CrTrHY(p, tr)).ToArray();
    LogLog(Pr, Cr);
    Tlabels.Add("Tr = " + tr);
}
Legend(Tlabels, Alignment.LowerLeft );
Xlabel("Pr"); Ylabel("CrTr");
Axis([0.1, 100, 0.01, 10]);
minorgrid = true;
SaveAs("Reduced-Compressibility-Hall-Yarborough-CCL-Math.png");

Blasius Boundary Layer and Howarth’s Transformation

Background

The Blasius Boundary Layer refers to a boundary layer of fluid in the vicinity of a flat plate that moves steadily in its own plane. This concept was first introduced by German mathematician Heinrich Blasius. This solution is important in the field of fluid dynamics, particularly in the area of laminar flow. In this solution, the flow velocity outside the boundary layer is assumed to be uniform. Inside the boundary layer, the fluid’s velocity changes from zero at the plate surface to the free stream velocity at the edge of the boundary layer. This concept plays a significant role in understanding and predicting the behavior of fluid flow in various engineering and scientific applications.

The solutions to the Blasius equation provide valuable insights into the behavior of fluid flow near a boundary. For instance, it demonstrates that the boundary layer thickness grows as the square root of the distance along the plate. Also, it reveals that the shear stress at the plate surface is proportional to the square root of the free stream velocity, among other observations.

The Blasius Boundary Layer solution, despite its simplifications, offers a good approximation for real-life engineering problems involving fluid flow over flat surfaces. This understanding is crucial in designing and optimizing various engineering systems, ranging from airfoils in aeronautics to heat exchangers in thermal power plants and thermal shield design on reusable rockets.

\[2f''' + f''f = 0\]

This equation can be solved by transforming it into a system of 1st order differential equations.

Let:

\[y_1 = f, y_2 = f', y_3 = f''\]

hence

\[\begin{split}\begin{array}{rcl} y'_1&=&f' = y_2, \\ y'_2&=&f'' = y_3, \\ y'_3&=&f''' = -0.5f''f = -0.5y_3y_1 \end{array}\end{split}\]

So, the system of first order differential equation is thus:

\[\begin{split}\begin{array}{rcl} y'_1&=&y_2, \\ y'_2&=&y_3, \\ y'_3&=&-0.5y_3y_1 \end{array}\end{split}\]

Subject to the following initial and terminal conditions.

\[y_1(0) = 0, y_2(0) = 0, y_2(\infty) = 1\]

To solve system of ordinary differential equation, you need the initial conditions, but when one of the initial conditions is missing, and we have a terminal condition instead, we can solve for the initial condition we do not have using the terminal condition we have, just like finding the root of a nonlinear function. But to evaluate the value of the function that you want to be zero, you have to perform the integration of the ode, using the guess given by the nonlinear solver and then return the function value which is the difference between the terminal value you obtained from the integration and the desired value. In this case the unknown value is \(y_3(0)\), and the function value is \(y_2(\infty) - 1\).

CCL-MathPythonMatlab

CCL-Math Implementation

// import libraries
using System;
using System.Collections.Generic;
using SepalSolver;
using static SepalSolver.Math;

// define variables
ColVec T = null; Matrix Y = null;

// define function
ColVec dydt(double t, ColVec y)
{
    double[] dy = [y[1], y[2], -0.5 * y[2] * y[0]];
    return dy;
}

// set time span
double[] tspan = [0, 6];

// define nonlinear function to shooting for terminal boundary
double fun(double y3_0)
{
    double[] y0 = [0, 0, y3_0];
    (T, Y) = Ode45(dydt, y0, tspan);
    return Y.LastRow[1] - 1;
}

// solve for unknown initial condition
Solvers.Result y3_0 = Solvers.FSolve(fun, 0.5);

// plot the result
Plot(T, Y, Linewidth: 2);
Legend( ["f", "f'", "f''"], Alignment.UpperLeft );
Axis([0, 6, 0, 2]); Xlabel("η"); Title("Blasius Boundary layer");
SaveAs("Blasius-Boundary-Layer-CCL-Math.png");

Howarth’s Transformation

Introducing the Howarth–Dorodnitsyn transformation (or Dorodnitsyn-Howarth transformation) is a density-weighted[4] coordinate transformation, which reduces variable-density flow conservation equations to simpler form (in most cases, to incompressible form), allows us to model the temperature within the boundary layer.

Introducing the self-similar variables

\[\eta = \sqrt{\frac{U}{2\nu_{\infty}x}} \int_0^y \frac{\rho}{\rho_{\infty}}dy, f(\eta) = \frac{\psi}{\sqrt{2\nu_{\infty}U_x}},~ \tilde{h}(\eta) = \frac{h}{h_{\infty}},~ \tilde{h}_w = \frac{h_w}{h_{\infty}},~ \tilde{\rho}_w = \frac{\rho}{\rho_{\infty}}\]

the equations reduce to:

\[2(\tilde{\rho}\tilde{\mu}f'')' + ff'' = 0,\]

\[(\tilde{\rho}\tilde{\mu}\tilde{h}')' + Pr f \tilde{h}' + Pr(\gamma - 1)M^2\tilde{\rho}\tilde{\mu}f''^2 = 0\]

where \(\gamma\) is the specific heat ratio and \(M = U/c_\infty\) is the Mach number[6], where \(c_\infty\) speed of sound[7]. The equation can be solved once \(\tilde{\rho} = \tilde{\rho}(\tilde{h})\) and \(\tilde{\mu} = \tilde{\mu}(\tilde{h})\) are specified.

The boundary conditions are:

\[f(0) = f’(0) = \theta(0) - \tilde{h}_w = f’(\infty)-1 = \tilde{h}(\infty)-1 = 0\]

The temperature inside the boundary layer will increase even though the plate temperature is maintained at the same temperature as ambient, due to dissipative heating and of course, these dissipation effects are only pronounced when the Mach number \(M\) is large.

Here we present the solution for the model when \(M = 0\) and \(5\).

CCL-MathPythonMatlab

CCL-Math Implementation

// import libraries
using System;
using System.Collections.Generic;
using SepalSolver;
using static SepalSolver.Math;

// define parameters
double rhomu_h, drhomu_h_eta, gamma, Pr, C;

// define functions and their derivatives
Func<double, double> rho, drhodh, mu, dmudh, rhomu;
Func<double, double, double> drhomu;

//define time span and intial guess
double[] tspan, y0, y35guess;

// define intexer for the unknwon initial conditions
Indexer I = new int[] { 1, 3 };

//define function for solution of howarth transformation
(ColVec, Matrix) HowarthTransform(double M)
{
    // assign parameters, functions anf their derivatives
    gamma = 1.4;
    Pr = 0.7;
    C = Pr * (gamma - 1) * M * M;
    rho = h => 1.0/h;
    drhodh = h => -1 / (h * h);
    mu = h => Pow(h, 2.0 / 3);
    dmudh = h => 2.0 / 3 * Pow(h, -1.0 / 3);
    rhomu = h => rho(h) * mu(h);
    drhomu = (h, dh) => (rho(h) * dmudh(h) + drhodh(h) * mu(h)) * dh;

    // define the differential equation
    ColVec dydt(double t, ColVec y)
    {
        rhomu_h = rhomu(y[3]);
        drhomu_h_eta = drhomu(y[3], y[4]);
        double[] dy = [y[1],
                       y[2],
                       -(2*drhomu_h_eta + y[0])*y[2]/(2*rhomu_h),
                       y[4],
                       -(drhomu_h_eta*y[4] + Pr*y[0]*y[4] + C*rhomu_h*y[2]*y[2])/rhomu_h ];
        return dy;
    }

    // set time span and intial guess
    tspan = [0, 5]; y35guess = [0.1, 0.2];
    ColVec T = null; Matrix Y = null;

    // define the nonlinear system to compute the initial condition
    ColVec fun(ColVec y35_0)
    {
        y0 = [0, 0, y35_0[0], 2, y35_0[1]];
        (T, Y) = Ode45(dydt, y0, tspan);
        return Y.LastRow[I].T - 1;
    }

    // solve for the unknown initial conditions
    Solvers.FSolve(fun, y35guess);
    return (T, Y);
}

// generator solution for M = 0 and plot
(ColVec T, Matrix Y) = HowarthTransform(0);
Plot(T, Y["", 1], "b", 2); hold = true;
Plot(T, Y["", 3] - 1, "r", 2);

// generator solution for M = 5 and plot
(T, Y) = HowarthTransform(5);
Plot(T, Y["", 1], "b", 2);
Plot(T, Y["", 3] - 1, "r", 2); hold = false;

// add legend, axis label and title
Legend( ["f'", "h-1"], Alignment.UpperRight);
Xlabel("η"); Title("Howarth Transformation");
Axis([0, 5, 0, 2]);
SaveAs("Howarth-Transformation-CCL-Math.png");

Pleiades System

The Pleiades, also known as the Seven Sisters (M45)[1], is a prominent open star cluster located in the constellation Taurus. It’s one of the closest and most easily visible star clusters to Earth[2], making it a favorite target for stargazers and a subject of fascination across cultures. The system of equations describing the motion of the stars in the cluster consists of 14 nonstiff second-order differential equations, which produce a system of 28 equations when rewritten in first-order form.

Celestial mechanics is basically an interplay between Newton’s law of gravitation \(F_i = \sum_{i \neq j} g \cfrac{m_i m_j}{||p_j - p_i||^2}d_{ij}\) and Newton’s second law of motion \(F_i = m_i\cfrac{d^2p_i}{dt^2}\).

The positions determine the gravitational forces acting on the bodies, but the net force on each of the bodies determines its acceleration (i.e. changes its position from the second order).

we examine this system in 2D, i.e. \(p_i = [x_i, y_i]\), \(d_{ij} = \cfrac{(p_j - p_i)}{r_{ij}}\) and \(r_{ij} = ||p_j - p_i||\)

The dynamics of the system can then be modelled as:

\[\cfrac{d^2p_i}{dt^2} = \sum_{i \neq j} g \cfrac{m_j(p_j - p_i)}{r_{ij}^3}\]

CCL-MathPythonMatlab

CCL-Math Implementation

// import libraries
using System;
using System.Collections.Generic;
using SepalSolver;
using static SepalSolver.Math;

// define masses
double[] m = [1, 2, 3, 4, 5, 6, 7];

// define function
ColVec pleiades(double t, ColVec q)
{
    double[] dqdt = new double[28];
    double x1, x2, y1, y2, dx, dy, r3;
    for (int i = 0; i < 7; i++)
    {
        // x- velocity of star i
        dqdt[i + 0] = q[i + 14];
        // y- velocity of star j
        dqdt[i + 7] = q[i + 21];
        x1 = q[i]; y1 = q[i + 7];
        for (int j = 0; j < 7; j++)
        {
            x2 = q[j]; y2 = q[j + 7];
            if (j != i)// The star does not attract itself
            {
                dx = x2 - x1; dy = y2 - y1;
                r3 = Pow(dx * dx + dy * dy, 1.5);
                //impact of star j on x-acceleration of star i
                dqdt[i + 14] += m[j] * dx / r3;
                //impact of star j on y-acceleration of star i
                dqdt[i + 21] += m[j] * dy / r3;
            }
        }
    }
    return dqdt;
}

double[] init = [3, 3,-1, -3, 2, -2, 2,
                 3, -3, 2, 0, 0, -4, 4,
                 0, 0, 0, 0, 0, 1.75, -1.5,
                 0, 0, 0, -1.25, 1, 0, 0];

Indexer I = new(0, 7), J = I + 7;
double[] tspan = Linspace(1, 15, 200);
var opts = Odeset( AbsTol: 1e-15, RelTol: 1e-13 );

(ColVec T, Matrix Y) = Ode89(pleiades, init, tspan, opts);
Plot(Y["", I], Y["", J], "--");
Title("Position of Pleiades Stars, Solved by ODE89");
Xlabel("X Position");
Ylabel("y Position");
SaveAs("Position-of-Pleiades-Stars-CCL-Math-Ode89.png");

we can add animation of the solution

CCL-MathPythonMatlab

CCL-Math Implementation

// import libraries
using System;
using System.Collections.Generic;
using SepalSolver;
using static SepalSolver.Math;


byte[] ImageGen(int i)
{
    Plot(Y["", I], Y["", J], "--"); hold = true;
    for (int j = 0; j < 7; j++) Scatter(Y[i, j], Y[i, j + 7], "fo", 20);
    hold = false; return GetImageBytes(700, 700);
}
AnimationMaker(ImageGen, "Position-of-Pleiades-Stars-CCL-Math-Ode89.gif", 10, 200);

One Transistor Amplifier

This example shows how to solve a stiff differential algebraic equation (DAE) that describes an electrical circuit. The one-transistor amplifier problem can be rewritten in semi-explicit form, but this example solves it in its original form \(Mu=\varphi(u)\) The problem includes a constant, singular mass matrix \(M\).

The transistor amplifier circuit contains six resistors, three capacitors, and a transistor.

• The initial voltage signal is \(U_e(t) = 0.4\sin(200\pi t)\).

• The operating voltage is \(U_b = 6\).

• The voltages at the nodes are given by \(U_i(t) (i = 1,2,3,4,5)\).

• The values of the resistors \(R_i(t) (i = 1,2,3,4,5)\). are constant, and the current through each resistor satisfies \(I = U/R\).

• The values of the capacitors \(C_i (i = 1,2,3)\) are constant, and the current through each capacitor satisfies \(I=C⋅dU/dt\).

The goal is to solve for the output voltage through node 5, \(U_5(t)\).

Using Kirchoff’s law to equalize the current through each node (1 through 5), you can obtain a system of five equations describing the circuit:

Node 1: \(C_1(U'_2 - U'_1) = (U_1 - U_e(t))/R_0\)

Node 2: \(C_1(U'_1 - U'_2) = (U_2 - U_b)/R_1 + U_2/R_1 + 0.01f(U_2 - U_3)\)

Node 3: \(-C_2U'_3 = U_3/R_3 - f(U_2 - U_3)\)

Node 4: \(C_3(U'_5 - U'_4) = (U_4 - U_b)/R_4 + 0.99f(U_2 - U_3)\)

Node 5: \(C_3(U'_4 - U'_5) = U_5/R_5\)

By extracting the coeeficients of the derivatives into a matrix, we have:

\[\begin{split}\begin{pmatrix} -c_{1} & c_{1} & 0 & 0 & 0 \\ c_{1} & -c_{1} & 0 & 0 & 0 \\ 0 & 0 & -c_{2} & 0 & 0 \\ 0 & 0 & 0 & -c_{3} & c_{3} \\ 0 & 0 & 0 & c_{3} & -c_{3} \end{pmatrix} \begin{pmatrix} U'_1 \\ U'_2 \\ U'_3 \\ U'_4 \\ U'_5 \end{pmatrix} = \begin{pmatrix} (U_1 - U_e(t))/R_0 \\ (U_2 - U_b)/R_1 + U_2/R_1 + 0.01f(U_2 - U_3) \\ U_3/R_3 - f(U_2 - U_3) \\ (U_4 - U_b)/R_4 + 0.99f(U_2 - U_3) \\ U_5/R_5 \end{pmatrix}\end{split}\]

CCL-MathPythonMatlab

CCL-Math Implementation

// import libraries
using System;
using System.Collections.Generic;
using SepalSolver;
using static SepalSolver.Math;

double pi = PI, Ub = 6, R0 = 1000, R15 = 9000, alpha = 0.99,
beta = 1e-6, Uf = 0.026, c1 = 1e-6, c2 = 2e-6, c3 = 3e-6;
Matrix Mass(double t, ColVec y) =>
   new double[,] { {-c1,  c1,  0,   0,   0 },
                   { c1, -c1,  0,   0,   0 },
                   { 0,   0,  -c2,  0,   0 },
                   { 0,   0,   0,  -c3,  c3},
                   { 0,   0,   0,   c3, -c3} };

ColVec dudt(double t, ColVec u)
{
    double Ue = 0.4 * Sin(200 * pi * t),
           f23 = beta * (Exp((u[1] - u[2]) / Uf) - 1);
    double[] du = [ -(Ue - u[0])/R0,
                  -(Ub/R15 - u[1]*2/R15 - (1-alpha)*f23),
                  -(f23 - u[2]/R15),
                  -((Ub - u[3])/R15 - alpha*f23),
                   u[4]/R15 ];
    return du;
}
double[] tspan = [0, 0.1];
double[] y0 = [0, Ub / 2, Ub / 2, Ub, 0];

Ode.Set options = new() { RelTol = 1e-3, MassType = Ode.MassType.Constant };

(ColVec T, Matrix Y) = Ode45a(dudt, Mass, y0, tspan, opts);
ColVec X = T, U5 = Y["", 4];
Scatter(X, 0.4 * Sin(200 * pi * X), "o"); hold = true;
Plot(X, U5, "--r"); hold = false;
Legend(["Input", "Output"], Alignment.UpperLeft);
Xlabel("Time t"); Ylabel("Solution y");
Title("One Transistor Amplifier DAE Problem-CCL-Math-DAE45");
SaveAs("One-Transistor-Amplifier-DAE-Problem-CCL-Math-DAE45.png");

Reference

1. Hairer, E., et al. Solving Ordinary Differential Equations I: Nonstiff Problems. 2nd rev. ed, Springer, 2009.

2. “Pleiades.” Wikipedia, 21 June 2021. Wikipedia, https://en.wikipedia.org/wiki/Pleiades.

3. Lateef Adewale Kareem (2025). Numerical Inversion of Laplace Transform (https://www.mathworks.com/matlabcentral/fileexchange/179769-numerical-inversion-of-laplace-transform), MATLAB Central File Exchange. Retrieved January 30, 2025.

4. https://en.wikipedia.org/wiki/Density

5. https://en.wikipedia.org/wiki/Specific_heat_ratio

6. https://en.wikipedia.org/wiki/Mach_number

7. https://en.wikipedia.org/wiki/Speed_of_sound