🎯要点

1. 根据温室模型，计算不同情景下辐射通量和评估能量平衡，构建复杂温室模型计算
2. 计算和绘图大气、海洋、陆地表面和海冰复合模型数据
3. 建立简单能量平衡情景模型，并根据模型计算释放温度和时滞，计算并绘制地面辐射和吸收光谱
4. 计算和模拟非散射辐射传输，求解史瓦西微分方程
   
   

🍇Python微分方程

如果因变量具有恒定的变化率：
$$\frac{dy}{dt}=C$$
其中 $C$ 是某个常数，您可以在 $f$ 函数中提供微分方程，然后使用以下代码使用此模型计算答案。该代码假设 0 和 10 之间有 100 个均匀间隔的时间， $y$ 的初始值为 6 ，变化率为 1.2 ：

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
from ode_helpers import state_plotterdef f(t, y, c):dydt = [c[0]]return dydttspan = np.linspace(0, 10, 100)
yinit = [6]
c = [1.2]sol = solve_ivp(lambda t, y: f(t, y, c), [tspan[0], tspan[-1]], yinit, t_eval=tspan, rtol = 1e-5)state_plotter(sol.t, sol.y, 1)

如果因变量的变化率是时间的函数，则可以轻松编码。例如，如果微分方程是某个二次函数，如下所示：
$$\frac{dy}{dt}=\alpha t^2+\beta t+\gamma$$
您可以使用此模型通过以下代码计算答案，假设 $y$ 的初始值 0 到 4 之间有 20 个均匀间隔的时间 𝑦 是 6，多项式由向量 [2, -6, 3] 定义：

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
from ode_helpers import state_plotterdef f(t, y, c):dydt = np.polyval(c, t)return dydttspan = np.linspace(0, 4, 20)
yinit = [6]
c = [2, -6, 3]sol = solve_ivp(lambda t, y: f(t, y, c), [tspan[0], tspan[-1]], yinit, t_eval=tspan, rtol = 1e-5)state_plotter(sol.t, sol.y, 1)

对于增长型，变化率取决于数量以及一些比例常数：
$$\frac{dy}{dt}=C\cdot y$$
其中 $C$ 又是某个常数。以下代码将计算上述模型的 3 秒内 25 个点的增长量，初始为 10，比例常数为 1.02：

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
from ode_helpers import state_plotterdef f(t, y, c):dydt = [c[0] * y[0]]return dydttspan = np.linspace(0, 3, 25)
yinit = [10]
c = [1.02]sol = solve_ivp(lambda t, y: f(t, y, c), [tspan[0], tspan[-1]], yinit, t_eval=tspan, rtol = 1e-5)state_plotter(sol.t, sol.y, 1)

通过确保微分函数返回每个变量的值，可以解决多变量系统。例如，在以下系统中，第一个变量的变化率仅取决于时间，而第二个变量则取决于时间和第一个变量：
$$\frac{d{y}_0}{dt}=\alpha \cos (\beta t)\frac{d{y}_1}{dt}=\gamma y_0+\delta t$$
该系统的微分函数 f 将有一个 2 元素列表作为输出。另外，如果您的系统具有多个因变量，请务必将初始条件放入列表中。例如，系统定义为：
$$\frac{d{y}_0}{dt}=4\cos (3t)\frac{d{y}_1}{dt}=-2y_0+0.5t$$
您可以使用以下脚本来求解 $y_0$ 和 $y_1$；代码假设 $y_0$ 从 0 开始，$y_1$ 从 -3 开始

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
from ode_helpers import state_plotterdef f(t, y, c):dydt = [c[0]*np.cos(c[1]*t), c[2]*y[0]+c[3]*t]return dydttspan = np.linspace(0, 5, 100)
yinit = [0, -3]
c = [4, 3, -2, 0.5]sol = solve_ivp(lambda t, y: f(t, y, c), [tspan[0], tspan[-1]], yinit, t_eval=tspan, rtol = 1e-5)state_plotter(sol.t, sol.y, 1)

系统必须仅以一阶微分方程的形式来写出。要求解具有高阶导数的系统，您首先要编写一个简单的一阶方程级联系统，然后在微分函数中使用它们。例如，假设您有一个以恒定抖动为特征的系统：
$$j=\frac{d^{3}y}{d{t}^3}=C$$
首先要做的是写出三个一阶微分方程来表示三阶方程：
$$\begin{array}{rl}& y_0=y\\ & \frac{d{y}_0}{dt}=\frac{dy}{dt}=y_1⟶\frac{d{y}_0}{dt}=y_1\\ & \frac{d{y}_1}{dt}=\frac{d^2{y}_0}{d{t}^2}=\frac{d^{2}y}{d{t}^2}=y_2⟶\frac{d{y}_1}{dt}=y_2\\ & \frac{d{y}_2}{dt}=\frac{d^2{y}_1}{d{t}^2}=\frac{d^3{y}_0}{d{t}^3}=\frac{d^{3}y}{d{t}^3}=j=C⟶\frac{d{y}_2}{dt}=C\end{array}$$
请注意导数如何级联，以便现在可以将恒定加加速度方程写为一组三个一阶方程。请注意，在该系统中，$y[0]$ 表示位置，$y[1]$ 表示速度，$y[2]$ 表示加速度。这种类型的级联系统在运动方程建模时经常出现。

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
from ode_helpers import state_plotterdef f(t, y, c):dydt = [y[1], y[2], c[0]]return dydttspan = np.linspace(0, 8, 50)
yinit = [6, 2, -4]
c = [1.3]sol = solve_ivp(lambda t, y: f(t, y, c), [tspan[0], tspan[-1]], yinit, t_eval=tspan, rtol = 1e-5)state_plotter(sol.t, sol.y, 1)

👉参阅、更新：计算思维 | 亚图跨际