我有一个修正的稳态一维传导方程：
a T（z）+ a z T“（z）+ B T”（z）= g T“（z）+ h
其中a、B、g、h > 0，z的范围为z1 = 0 m且z2 > 0 m，T的范围为T1和T2（其中T2 > T1 > 0 K）。
注意，a、B、g和h是与一个或多个其它物理常数相关的物理常数（即，不是z或t的函数），并且z1、z2、T1、T2也是常数。
一般来说，方程的左侧也可以有一个**j T '（t）**项（其中j > 0是一个常数），尽管在这个物理系统的上下文中，在相对较短的时间尺度内达到稳态，所以我们可以假设dT/dt = 0 K/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数值，蓝色实线=分析）：
我想用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 + a z）/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 '（s）z2^2 + B T'（s）z2 = g T“（s）+ h z2^2
下面是相关的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方法解决此问题，并使用子方法Newton。
我在网上寻找solve_bvp（）的替代Python函数/方法，特别是接近Mathematica的NDSolve[]或LSODA，还没有找到合适的。
据我所知，Mathematica的NDSolve[]方法LSODA是从Fortran或C算法修改而来的。
我发现Python确实有一个LSODA函数，它使用Fortran算法。
然而，它的设置方式（与Python的solve_ivp（）和odeint（）函数沿着）需要一个T（s）配置文件的初始数组（这正是我试图解决的问题），并且随着时间的推移而演变（而我的情况是与时间无关的）。
有没有一种相对快速的方法可以使用Python的scipy.integrate.LSODA（）或scipy.integrate.solve_ivp（method = 'LSODA'）或scipy.integrate.odeint（）或其他Python方法/函数来解决手头的BVP问题（并返回与Mathematica的NDSolve[]相同的解决方案）？
将**j T '（t）**项引入到上面给出的第一个方程的左手边，并可能将t无量纲化为物理系统的适当尺度（例如，tau = t/tM，其中tM的数量级为~0.1-1百万年），是否允许我通过使用温度分布T（s）的初始猜测（不是实际温度分布T（s））来求解BVP？
如果需要测试，我可以提供分析解决方案，和/或，Python和Mathematica的s和T输出的数组，用于上面的特定参数/常量。
请注意，我已经尝试过通过无量纲化T来求解这个系统（chi = [T-T1]/[T2 - T1]，边界条件chi（s = 0）= 0和chi（s = 1）= 1），以及z和T都没有无量纲化。这并没有改变我的Python solve_bvp（）温度配置文件解决方案。
2条答案
按热度按时间3ks5zfa01#
你要找的是一种“射击方法”。从本质上讲，您将BVP转换为IVP，并尝试找到一阶导数的初始值，以便在区间结束时解决方案满足边界条件。Here is a python walk-through.你可以使用任何你想要的好的IVP求解器。
你应该知道，射击方法可能不工作，在所有的一些问题。它似乎与你的工作，但它是你应该记住的东西。
kr98yfug2#
在看到@Mario S之前。Mommer的建议，我偶然发现的Finite Difference Method。
对于热方程：
ρ cp T“（t）= k T”（z）
其中alpha = k/（rho cp），T（z1 = 0）= T1，T（z2）= T2，我们可以在Python中通过调整有限差分方法的online example来解决它，使用以下代码：
这一做法产生了以下结果：
有限差分法可能适用于我修改后的热方程，使用这些中心差分近似：
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
然而，马里奥的射击方法的建议显着提高了准确性（与上述有限差分方法的实现相比）。
对于该稳态热方程：
k T“（z）+ Q = 0
对于T（z1 = 0）= T1和T（z2）= T2边界条件，我们可以通过修改this Shooting Method example来解决Python中的问题，使用以下代码：
这一做法产生了以下结果：
将这个Shooting Method Python实现适应我修改过的传导方程：
这将产生以下结果：