代码之家 › 专栏 › 技术社区 › Canada709

如何使用类似Mathematica的NDSolve[](方法LSODA和子方法Newton)的函数在Python中求解边值问题?

numerical-methods wolfram-mathematica scipy python

-1

Canada709 · 技术社区 · 1 年前

我有一个修正的稳态一维传导方程:

a T(z)+a z T’(z)+b T’(z)=g T’(兹)+h

其中a,b,g,h>z的范围从z1=0m到z2>并且T的范围为T1和T2(其中T2>T1>0K)。

注意,a、b、g和h是与一个或多个其他物理常数有关的物理常数(即,不是z或t的函数),并且z1、z2、T1T2也是常数。

一般来说,该方程也可以有一个 j T’(T) 左边的项(其中j>0是常数),然而,在这个物理系统的上下文中,在相对较短的时间尺度内达到稳态,因此,我们可以假设dT/dT=0K/s。

我已经使用Python和Mathematica数值求解了这个边值问题(BVP),对于与这个物理系统相关的各种常数的一些特定值,并且边界条件为T(z=z1)=T1和T(z=z2)=T2。

为了解决一个不同但相关的问题,我还对z参数进行了无量纲化,s=z/z2。

然而,绘制Python的scipy.integrate.solve_bvp()和Mathematica的NDSolve[]的解决方案显示了在s=0到s=1(对应于z=0到z=z2)范围内的不同温度分布T(s)。

以下是一个这样的示例图(红点=Python数值,绿点=Mathematica数值,蓝实线=Analytical):

我想用Python来解决这个问题,因为我使用这个温度配置文件的超长代码的其余部分都是用Python编写的。我只是出于测试和调试的目的在Mathematica上尝试了一下。

通过求解该通解的积分常数c1和c2,存在该问题的解析精确解:

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+az)/Sqrt[2 g a]]

然而,一般来说,这对我的问题来说不是一个可行的解决方案,因为在这个物理系统的一些相关参数空间(z1,z2,T1,T2,a,b,g,&h)中,解析解由于几个不同的原因而失败。

上面显示的示例图中,分析解决方案不会失败。

该图还显示,Mathematica的NDSolve[]解决方案与分析解决方案匹配,而Python的solve_bvp()则不匹配。

Python&Mathematica成功地拟合了T(z=z1)=T1和T(z=z2)=T2的边界条件。

这是z非维化方程,使用s=z/z2,我用Python和Mathematica对其进行了数值求解:

a T(s)z2^2+a s T’

这是相关的Python代码:

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()

这是相关的Mathematica代码:

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[]"]

NDSolve[]的Mathematica Trace[]表明,NDSolve[]使用了LSODA方法来解决这个问题,并使用了子方法牛顿。

我在网上寻找solve_bvp()的替代Python函数/方法,特别是接近Mathematica的NDSolve[]或LSODA的函数/方法。

据我所知,Mathematica的NDSolve[]方法LSODA是从同名的Fortran或C算法中修改而来的。

我发现Python确实有一个LSODA函数,它使用Fortran算法。

然而,它的设置方式(以及Python的solve_ivp()和odeint()函数)需要一个T(s)配置文件的初始数组(这正是我试图解决的问题),并随着时间的推移而不断发展(而我的情况与时间无关)。

有没有一种相对快速的方法可以使用Python的scipy.integer.LSODA()、scipy.inintegre.solve_ivp(method='LSODA')、scipy.integer.odeint()或其他Python方法/函数来解决手头的BVP问题(并返回与Mathematica的NDSolve[]相同的解决方案)?

会介绍 j T’(T) 上面给出的第一个方程左手边的项,以及可能将t非尺寸化为物理系统的适当标度(例如,tau=t/tM,其中tM的数量级约为10-10万年),允许我通过使用不是实际温度分布t(s)的温度分布t的初始猜测来求解BVP?

如果需要测试和/或Python和amp;上面特定参数/常数的Mathematica。

注意,我已经尝试过通过对T进行非维化来求解这个系统(其中chi=[T-T1]/[T2-T1],边界条件为chi(s=0)=0和chi(s=1)=1),并且z和T都没有进行非维。这并没有改变我的Python solve_bvp()温度配置文件解决方案。

0 回复 | 直到 1 年前

Mario S. Mommer 1 年前

你要找的是 射击方法 本质上,您将BVP转换为IVP,并尝试找到一阶导数的初始值,以便在区间结束时解满足边界条件。 Here is a python walk-through. 您可以使用任何您想要的好的IVP解算器。

您应该意识到,在出现某些问题时,拍摄方法可能根本不起作用。它似乎和你的一样有效,但这是你应该记住的。

Canada709 1 年前

在看到@Mario S.Mommer的建议之前,我遇到了 Finite Difference Method 。

对于热量方程:

ρcp T'(T)=k T'(z)

当alpha=k/(rho-cp),T(z1=0)=T1,T(z 2)=T2时,我们可以在Python中通过调整 online example 对于有限差分法,使用以下代码:

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()

这种方法产生以下结果:

有限差分法可能适用于我修改的热方程,使用这些中心差分近似值:

T'(z_i)~T(z_(i+1))-T(z_

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

尽管如此,马里奥对射击方法的建议显著提高了精度(与上述有限差分法的实现相比)。

对于这个稳态热方程:

k T''(z)+Q=0

当T(z1=0)=T1和T(z2)=T2边界条件时,我们可以通过调整 this Shooting Method example ,使用此代码:

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()

这种方法产生以下结果:

将此Shooting Method Python实现应用于我修改的传导方程:

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()

这将产生以下内容: