🥭Python气象辐射光谱能量平衡模型

Python气象辐射光谱能量平衡模型 | 亚图跨际亚图跨际 on Notion

🏈指点迷津 | Brief

🎯要点

1. 根据温室模型，计算不同情景下辐射通量和评估能量平衡，构建复杂温室模型计算

2. 计算和绘图大气、海洋、陆地表面和海冰复合模型数据

3. 建立简单能量平衡情景模型，并根据模型计算释放温度和时滞，计算并绘制地面辐射和吸收光谱

4. 计算和模拟非散射辐射传输，求解史瓦西微分方程

🍪语言内容分比

🍇Python微分方程

如果因变量具有恒定的变化率：

$$\frac{d y}{d t}=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_plotter
 ​
 def f(t, y, c):
     dydt = [c[0]]
     return dydt
 ​
 tspan = 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{d y}{d t}=\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_plotter
 ​
 def f(t, y, c):
     dydt = np.polyval(c, t)
     return dydt
     
 tspan = 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{d y}{d t}=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_plotter
 ​
 def f(t, y, c):
     dydt = [c[0] * y[0]]
     return dydt
 ​
 tspan = 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}{d t}=\alpha \cos (\beta t) \quad \frac{d y_1}{d t}=\gamma y_0+\delta t$$

该系统的微分函数 f 将有一个 2 元素列表作为输出。另外，如果您的系统具有多个因变量，请务必将初始条件放入列表中。例如，系统定义为：

$$\frac{d y_0}{d t}=4 \cos (3 t) \quad \frac{d y_1}{d t}=-2 y_0+0.5 t$$

您可以使用以下脚本来求解 $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_plotter
 ​
 def f(t, y, c):
     dydt = [c[0]*np.cos(c[1]*t), c[2]*y[0]+c[3]*t]
     return dydt
 ​
 tspan = 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{aligned} & y_0=y \\ & \frac{d y_0}{d t}=\frac{d y}{d t}=y_1 \quad \longrightarrow \quad \frac{d y_0}{d t}=y_1 \\ & \frac{d y_1}{d t}=\frac{d^2 y_0}{d t^2}=\frac{d^2 y}{d t^2}=y_2 \quad \longrightarrow \quad \frac{d y_1}{d t}=y_2 \\ & \frac{d y_2}{d t}=\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 \quad \longrightarrow \quad \frac{d y_2}{d t}=C \end{aligned}$$

请注意导数如何级联，以便现在可以将恒定加加速度方程写为一组三个一阶方程。请注意，在该系统中，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_plotter
 ​
 def f(t, y, c):
     dydt = [y[1], y[2], c[0]]
     return dydt
 ​
 tspan = 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)