How To Solve A Boundary Value Problem In Python Using A Function Like Mathematica's NDSolve[] (method LSODA and sub-method Newton)?

I have a modified steady-state 1-D conduction equation:

a T(z) + a z T'(z) + b T'(z) = g T''(z) + h

where a, b, g, h > 0, z ranges from z1 = 0 m and z2 > 0 m, and T ranges from T1 and T2 (with T2 > T1 > 0K).

Note that a, b, g, and h are physical constants (i.e., not functions of z nor t), relating to one or more other physical constants, and, z1, z2, T1 T2 are constants as well.

In general, the equation could also have a j T'(t) term on the left-hand-side (where j > 0 is a constant), though, in the context of this physical system, a steady-state is reached in a relatively short timescale, so, we can assume dT/dt = 0 K/s.

I have numerically solved this boundary value problem (BVP) using both Python and Mathematica, for some specific values of the various constants relevant for this physical system, and, with boundary conditions of T(z = z1) = T1 and T(z = z2) = T2.

I also non-dimensionalized the z-parameter, with s = z/z2, for sake of solving a different but related problem.

However, plotting the solutions from Python's scipy.integrate.solve_bvp() and Mathematica's NDSolve[] shows different temperature profiles T(s) for the range of s = 0 to s = 1 (corresponding to z = 0 to z = z2).

Here is one such example plot (red dotted = Python numerical, green dashed = Mathematica numerical, blue solid = Analytical):

I want to solve the problem with Python, because the rest of my very-long code which uses this temperature profile is written in Python. I only tried it on Mathematica for testing and debugging sake.

There exists an analytical-exact solution to this problem, by solving for the integration constants c1 and c2 of this general solution:

T(z) = h/a + c2 exp[z[2b + az]/2g] + Sqrt[pi/(2 g a^3)] (h b + g a c1) Exp[(b + a z)^2/(2 g a)] Erf[(b + a z)/Sqrt[2 g a]]

However, in general, this isn't a viable solution for my problem, because the analytical solution fails for a few different reasons in some of the relevant parameter spaces (of z1,z2,T1,T2,a,b,g,&h) for this physical system.

The example plot shown above is one in which the analytical solution does not fail.

The plot also shows that Mathematica's NDSolve[] solution matches the analytical solution, while Python's solve_bvp() does not.

Both Python & Mathematica succeed in fitting the boundary conditions of T(z = z1) = T1 and T(z = z2) = T2.

This is the z-non-dimensionalized equation, using s = z/z2, which I numerically solved with both Python and Mathematica:

a T(s) z2^2 + a s T'(s) z2^2 + b T'(s) z2 = g T''(s) + h z2^2

This is the relevant Python code:

import numpy as np
from scipy.integrate import solve_bvp
import matplotlib.pyplot as plt

def func(x, y, z2, a, b, g, h):

    second_deriv_arr = np.array([])
    for i in np.arange(len(x)):
        s_val = x[i]
        T_val = y[0][i]   
        dT_ds_val = -y[1][i]

        second_deriv_val = -1.*(h/g)*(z2**2) + (a/g)*T_val*(z2**2) + (a/g)*s_val*dT_ds_val*(z2**2) + (b/g)*dT_ds_val*z2
        second_deriv_arr = np.append(second_deriv_arr, second_deriv_val)

    return np.vstack((y[1], second_deriv_arr))


def bc(ya, yb, T1, T2):

    return np.array([ya[0] - T2, yb[0] - T1])


def get_res_func(z2, T1, T2, a, b, g, h, tolerance = 1e-4, step_size_km = 0.5):

    #z1 is 0 m for this physical system
    #z2 is ~10-100km for this physical system
    step_size = step_size_km*1000.
    gap1 = int(z2/step_size) + 1 
    x_in = np.linspace(0, 1, gap1)
    x_length = len(x_in)

    y_in = np.zeros((2, x_length))
    y_in[0] = np.linspace(T2, T1, x_length)
    y_in[1] = -(T2 - T1)/1 #would be division by z2 here if is the system wasn't non-dimensionalized

    res_func = solve_bvp(lambda x,y: func(x, y, z2, a, b, g, h), lambda ya,yb: bc(ya, yb, T1, T2), x_in, y_in, tol = tolerance)
    #could check res_func.success Truth value here

    s_arr = x_in #could multiply by z2 to re-deminsionalize
    T_arr = np.flip(res_func.sol(x_in)[0])

    return s_arr, T_arr, res_func

#constant values close to the ones used for the above plot
z1 = 0.
z2 = 20.*1000
T1 = 100.
T2 = 1500.
a = 7/1.5768e8
b = 7./7884
g = 3.
h = 7./140160

s_arr, T_arr, res_func = get_res_func(z2, T1, T2, a, b, g, h)

plt.plot(s_arr, T_arr, 'o', color = 'r', alpha = 0.95, markersize = 6, markeredgecolor = 'none', label = 'Python solve_bvp()')
plt.xlim(xmin = 0., xmax = 1.)
plt.ylim(ymin = T1, ymax = T2)
plt.xlabel("s", fontsize = 16)
plt.ylabel("Temperature [K]", fontsize = 16)
plt.legend(fontsize = 12)
plt.show()

This is the relevant Mathematica code:

LHS[z2_,a_,b_,g_,h_]:=a*z2*z2*T[s]+a*s*z2*z2*T'[s]+b*z2*T'[s]
RHS[z2_,a_,b_,g_,h_]:=g*T''[s]+h*z2*z2

#constant values close to the ones used for the above plot
z1 = 0.
z2 = 20.*1000
T1 = 100.
T2 = 1500.
a = 7/1.5768e8
b = 7./7884
g = 3.
h = 7./140160

sol = NDSolve[{LHS[z2,a,b,g,h] == RHS[z2,a,b,g,h], T[0]==T1, T[1]==T2}, T, {s, 0, 1}]
(*Trace[NDSolve[{LHS[z2,a,b,g,h] == RHS[z2,a,b,g,h], T[0]==T1, T[1]==T2}, T, {s, 0, 1}],_[NDSolve`MethodData[___]],TraceInternal->True]*)

Plot[[T[k] /. sol, {k, 0, 1},PlotRange -> {{0, 1}, {100, 1500}}, PlotStyle -> {Dashed, Green, Thick}, AxesLabel -> {"s", "T [K]"}, PlotLegends -> "-- NDSolve[]"]

The Mathematica Trace[] of NDSolve[] shows that NDSolve[] used the method LSODA for this problem, with sub-method Newton.

I looked online for alternative Python functions/methods for solve_bvp(), specifically to be close to Mathematica's NDSolve[] or LSODA, and haven't yet come across one that seems appropriate.

From what I have read, Mathematica's NDSolve[] method LSODA is modified from either a Fortran or C algorithm of the same name.

I have found that Python does have a LSODA function, which uses the Fortran algorithm.

However, the way it is set up (along with Python's solve_ivp() and odeint() functions) requires an initial array of the T(s) profile (which is precisely what I'm trying to solve for), and, evolves it over time (whereas my case is time-independent).

Is there a relatively quick way to use Python's scipy.integrate.LSODA(), or scipy.integrate.solve_ivp(method = 'LSODA'), or scipy.integrate.odeint(), or some other Python method/function to solve the BVP problem at hand (and return the same solution as Mathematica's NDSolve[])?

Would introducing the j T'(t) term to the left-hand-side of the first equation given above, and possibly non-dimensionilizing t to an appropriate scaling for the physical system (e.g., tau = t/tM, where tM is of order ~0.1-1 million years), allow for me to solve the BVP, by using an initial guess of the temperature profile T(s) that is not the actual temperature profile T(s)?

I can provide the analytical solution if needed for testing, and/or, arrays of the s and T outputs of both Python & Mathematica for the particular parameters/constants above.

Note that I have already tried solving this system with by non-dimensionilizing T as well (with chi = [T - T1]/[T2 - T1], and boundary conditions of chi(s = 0) = 0 and chi(s = 1) = 1), as well as with neither z nor T non-dimensionlized. This did not change my Python solve_bvp() Temperature profile solution.

Before seeing @Mario S. Mommer 's suggestion, I came across the Finite Difference Method.

For the heat equation:

rho cp T'(t) = k T''(z)

with alpha = k/(rho cp), T(z1 = 0) = T1, T(z2) = T2, we can solve it in Python by adapting an online example of the finite differences method, using this code:

import numpy as np
import matplotlib.pyplot as plt


#assumes that z1 = 0
def find_T_arr(T1, T2, z2, alpha, t_final, dt, n = int(20)):

    dz = z2/n
    z_arr = np.linspace(dz, z2 - dz, n)

    T_guess = (T2 - T1)/2
    T_arr = np.ones(n)*T_guess
    dTdt_arr = np.empty(n)

    t_arr = np.arange(0, t_final + dt, dt) 

    for j in range(1, len(t_arr)):
        for i in range(1, n-1):
            dTdt_arr[i] = alpha*((T_arr[i+1]-T_arr[i])/dz**2 -(T_arr[i] - T_arr[i-1])/dz**2)
        dTdt_arr[0] = alpha*((T_arr[1]-T_arr[0])/dz**2 -(T_arr[0] - T1)/dz**2)
        dTdt_arr[n-1] = alpha*((T2-T_arr[n-1])/dz**2 -(T_arr[n-1] - T_arr[n-2])/dz**2)
        T_arr = T_arr + dTdt_arr*dt

    z_arr = np.insert(z_arr, 0, 0.) #assumes z1 = 0
    z_arr = np.append(z_arr, z2)
    T_arr = np.insert(T_arr, 0, T1)
    T_arr = np.append(T_arr, T2)

    return z_arr, T_arr


#assumes z1 = 0 and at steady-state (e.g., dT/dt = 0)
def get_T_exact(z_arr, T1, T2, z2):

    Texact_arr = np.array([])
    for z_val in z_arr:

        T_val = T1 + ((T2 - T1)/(z2 - 0.))*z_val
        Texact_arr = np.append(Texact_arr, T_val)

    return Texact_arr


#z1 = 0. #assumed in the functions
z2 = 20.*1000
T1 = 100.
T2 = 1500.
alpha = 3/(1000*2800)
dt = 0.005*(365*24*60*60)*1e6
t_final = 10*(365*24*60*60)*1e6 #could be smaller

z_arr, T_arr = find_T_arr(T1, T2, z2, alpha, t_final, dt, n = int(30))
Texact_arr = get_T_exact(z_arr, T1, T2, z2)

fig, [ax1, ax2] = plt.subplots(nrows = 1, ncols = 2, sharex = False, sharey = False, figsize = (12, 10))

ax1.plot(z_arr/1000., Texact_arr, color = 'blue', linewidth = 2, label = 'Exact')
ax1.plot(z_arr/1000., T_arr, 'o', color = 'red', alpha = 0.95, markersize = 6, markeredgecolor = 'none', label = 'Finite Differences')
ax1.set_xlim(xmin = 0., xmax = z2/1000.)
ax1.set_ylim(ymin = 0., ymax = T1 + T2)
ax1.set_xlabel("Depth [km]", fontsize = 16)
ax1.set_ylabel("Temperature [K]", fontsize = 16)
ax1.hlines(y = T1, xmin = 0., xmax = z2/1000., linestyle = '--', color = 'grey', linewidth = 1)
ax1.hlines(y = T2, xmin = 0., xmax = z2/1000., linestyle = '--', color = 'grey', linewidth = 1)
ax1.legend(fontsize = 12, framealpha = 0.95)

ax2.plot(z_arr[1:-1]/1000., Texact_arr[1:-1] - T_arr[1:-1], color = 'purple', linewidth = 2)
ax2.set_xlim(xmin = 0., xmax = z2/1000.)
ax2.set_ylim(ymin = -2., ymax = 2.)
ax2.set_xlabel("Depth [km]", fontsize = 16)
ax2.set_ylabel("Temperature Residuals [K]", fontsize = 16)
ax2.hlines(y = 0., xmin = 0., xmax = z2/1000., linestyle = '--', color = 'grey', linewidth = 1)

fig.show()

That approach produces the following:

The Finite Differences method could probably be adapted to my modified heat equation, using these Central Differences approximations:

T'(z_i) ~ T(z_(i+1)) - T(z_(i-1))/(2*dz)

T''(z_i) ~ (T(z_(i+1)) - T(z_i))/dz^2 - (T(z_(i)) - T(z_(i-1)))/dz^2

Though, Mario's suggestion of the Shooting Method increases the accuracy significantly (in comparison to the above Finite Differences Method implementation).

For this steady-state heat equation:

k T''(z) + Q = 0

with T(z1 = 0) = T1 and T(z2) = T2 boundary conditions, we can solve the problem in Python by adapting this Shooting Method example, using this code:

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
from scipy.optimize import fsolve


def F_func(z, s, Q, k):
    F_val = [s[1], -Q/k]
    #z is the depth
    #s[0] is T(z)
    #s[1] is beta(z) = dT(z)/dz
    #F_val[1] is T''(z)
    return F_val


#assumes z1 = 0
def find_actual_beta0(beta0_val, Q, k, T1, T2, z2, z_eval):

    sol = solve_ivp(lambda z,s: F_func(z, s, Q, k), t_span = [0., z2], y0 = [T1, beta0_val], t_eval = z_eval)
    T_arr = sol.y[0]

    return T_arr[-1] - T2


def find_T_arr(T1, T2, z2, Q, k, n = int(10)):

    z_eval = np.linspace(0, z2, n) #assumes z1 = 0
    beta0_guess = (T2 - T1/2.)/(z2 - 0.) #assumes z1 = 0

    beta0 = fsolve(lambda beta0_val: find_actual_beta0(beta0_val, Q, k, T1, T2, z2, z_eval), x0 = beta0_guess)
    sol = solve_ivp(lambda z,s: F_func(z, s, Q, k), t_span = [0, z2], y0 = [T1, beta0], t_eval = z_eval)

    z_arr = sol.t
    T_arr = sol.y[0]

    return z_arr, T_arr


def get_Texact_val(z, z2, k, Q, c1, c2):

    #assumes z1 = 0, and 0 <= z <= z2
    Texact_val = -(Q*(z**2))/(2*k) + c1 + c2*z

    return Texact_val


def get_Texact_arr(z_arr, z2, T1, T2, k, Q):

    c1 = T1 #assumes z1 = 0 
    c2 = (z2*Q)/(2*k) + (T2 - T1)/z2 #assumes z1 = 0

    Texact_arr = np.array([])
    for z_val in z_arr:

        T_val = get_Texact_val(z_val, z2, k, Q, c1, c2)
        Texact_arr = np.append(Texact_arr, T_val)

    return Texact_arr


#z1 = 0. #assumed in the functions
z2 = 20e3
T1 = 100.
T2 = 1500.
Q = 0. #possible values for the system are 0 <= Q <= 1e-5 W/m^3
k = 3.

BC_z = np.array([0., z2])
BC_T = np.array([T1, T2])

z_arr, T_arr = find_T_arr(T1, T2, z2, Q, k)
Texact_arr = get_Texact_arr(z_arr, z2, T1, T2, k, Q)

fig, [ax1, ax2] = plt.subplots(nrows = 1, ncols = 2, sharex = False, sharey = False, figsize = (12, 10))

ax1.plot(z_arr/1000., Texact_arr, color = 'blue', linewidth = 2, label = 'Exact')
ax1.plot(BC_z/1000., BC_T, 'o', color = 'yellow', alpha = 0.95, markersize = 12, markeredgecolor = 'none', label = 'Boundary Conditions')
ax1.plot(z_arr/1000., T_arr, 'o', color = 'r', alpha = 0.95, markersize = 6, markeredgecolor = 'none', label = 'Shooting Method')
ax1.set_xlabel("Depth [km]", fontsize = 16)
ax1.set_ylabel("Temperature [K]", fontsize = 16)
ax1.set_xlim(xmin = 0., xmax = z2/1000.)
ax1.set_ylim(ymin = 0., ymax = T1 + T2)
ax1.legend(fontsize = 12, framealpha = 0.95)
ax1.hlines(y = T1, xmin = 0., xmax = z2/1000., linestyle = '--', color = 'grey', linewidth = 1)
ax1.hlines(y = T2, xmin = 0., xmax = z2/1000., linestyle = '--', color = 'grey', linewidth = 1)

ax2.plot(z_arr/1000., Texact_arr - T_arr, color = 'purple', linewidth = 2)
ax2.set_xlim(xmin = 0., xmax = z2/1000.)
ax2.set_ylim(ymin = 0., ymax = 2.5e-12)
ax2.set_xlabel("Depth [km]", fontsize = 16)
ax2.set_ylabel("Temperature Residuals [K]", fontsize = 16)

fig.show()

That approach produces the following:

Adapting this Shooting Method Python implementation to my modified conduction equation:

import numpy as np
from math import erf as err_func
import math
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
from scipy.optimize import fsolve


def F_func(z, s, a, b, g, h):
    F_val = [s[1], -h/g + (a/g)*s[0] + (a*z/g)*s[1] + (b/g)*s[1]]
    #z is the depth
    #s[0] is T(z)
    #s[1] is beta(z) = dT(z)/dz
    #F_val[1] is T''(z)
    return F_val


#assumes z1 = 0
def find_actual_beta0(beta0_val, a, b, g, h, T1, T2, z2, z_eval):

    sol = solve_ivp(lambda z,s: F_func(z, s, a, b, g, h), t_span = [0., z2], y0 = [T1, beta0_val], t_eval = z_eval, method = 'Radau')
    #tried method = 'RK45' and 'LSODA', found that 'Radau' had better accuracy for this problem
    T_arr = sol.y[0]

    return T_arr[-1] - T2


def find_T_arr(T1, T2, z2, a, b, g, h, n = int(10)):

    z_eval = np.linspace(0, z2, n) #assumes z1 = 0
    beta0_guess = (T2 - T1/2.)/(z2 - 0.) #assumes z1 = 0

    beta0 = fsolve(lambda beta0_val: find_actual_beta0(beta0_val, a, b, g, h, T1, T2, z2, z_eval), x0 = beta0_guess)
    sol = solve_ivp(lambda z,s: F_func(z, s, a, b, g, h), t_span = [0, z2], y0 = [T1, beta0], t_eval = z_eval, method = 'Radau')
    #tried method = 'RK45' and 'LSODA', found that 'Radau' had better accuracy for this problem

    z_arr = sol.t
    T_arr = sol.y[0]

    return z_arr, T_arr


def get_Texact_val(z, c1, c2, a, b, g, h):

    #assumes z1 = 0, and 0 <= z <= z2
    Texact_val = c2*np.exp((z*(2*b + a*z))/(2*g)) + h/a + (math.sqrt(math.pi/(2*g*(a**3.))))*(h*b + g*a*c1)*err_func((b + a*z)/(math.sqrt(2*g*a)))*np.exp(((b + a*z)**2)/(2*g*a))

    return Texact_val


def get_Texact_arr(z_arr, z2, T1, T2, a, b, g, h):

    exp1 = np.exp(((2.*b + a*z2)/(2*g))*z2)
    exp2 = np.exp(((b + a*z2)**2)/(2.*g*a))
    exp3 = np.exp((b**2)/(2.*g*a))
    erf1 = err_func((b + a*z2)/(math.sqrt(2.*g*a)))
    erf2 = err_func(b/(math.sqrt(2.*g*a)))

    c1 = (a*b*((-erf1)*exp2 + erf2*exp1*exp3)*h + (-1. + exp1)*math.sqrt(a**3*g)*h*math.sqrt(2/math.pi) + a*math.sqrt(a**3*g)*math.sqrt(2/math.pi)*(T2 - exp1*T1))/(a**2*(erf1*exp2 - erf2*exp1*exp3)*g) 
    c2 = (erf2*exp3*(h - a*T2) + erf1*exp2*(-h + a*T1))/(a*(erf1*exp2 - erf2*exp1*exp3))

    Texact_arr = np.array([])
    for z_val in z_arr:

        T_val = get_Texact_val(z_val, c1, c2, a, b, g, h)
        Texact_arr = np.append(Texact_arr, T_val)

    return Texact_arr


#z1 = 0. #assumed in the functions
z2 = 20e3
T1 = 100.
T2 = 1500.
a = 7/1.5768e8
b = 7./7884
g = 3.
h = 7./140160

BC_z = np.array([0., z2])
BC_T = np.array([T1, T2])

z_arr, T_arr = find_T_arr(T1, T2, z2, a, b, g, h, n = int(20))
Texact_arr = get_Texact_arr(z_arr, z2, T1, T2, a, b, g, h)

fig, [ax1, ax2] = plt.subplots(nrows = 1, ncols = 2, sharex = False, sharey = False, figsize = (12, 10))
fig.subplots_adjust(wspace = 0.3)

ax1.plot(z_arr/1000., Texact_arr, color = 'blue', linewidth = 2, label = 'Exact')
ax1.plot(BC_z/1000., BC_T, 'o', color = 'yellow', alpha = 0.95, markersize = 12, markeredgecolor = 'none', label = 'Boundary Conditions')
ax1.plot(z_arr/1000., T_arr, 'o', color = 'r', alpha = 0.95, markersize = 6, markeredgecolor = 'none', label = 'Shooting Method')
ax1.set_xlabel("Depth [km]", fontsize = 16)
ax1.set_ylabel("Temperature [K]", fontsize = 16)
ax1.set_xlim(xmin = 0., xmax = z2/1000.)
ax1.set_ylim(ymin = 0., ymax = T1 + T2)
ax1.legend(fontsize = 12, framealpha = 0.95)
ax1.hlines(y = T1, xmin = 0., xmax = z2/1000., linestyle = '--', color = 'grey', linewidth = 1)
ax1.hlines(y = T2, xmin = 0., xmax = z2/1000., linestyle = '--', color = 'grey', linewidth = 1)

ax2.plot(z_arr/1000., Texact_arr - T_arr, color = 'purple', linewidth = 2)
ax2.set_xlim(xmin = 0., xmax = z2/1000.)
ax2.set_ylim(ymin = -0.0125, ymax = 0.0075)
ax2.set_xlabel("Depth [km]", fontsize = 16)
ax2.set_ylabel("Temperature Residuals [K]", fontsize = 16)
ax2.hlines(y = 0., xmin = 0., xmax = z2/1000., linestyle = '--', color = 'grey', linewidth = 1)

fig.show()

That produces the following: