如何绘制与现有下降曲线反向增长的第二条Python曲线

原代码及运行结果 #

原代码 #

import scipy
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.collections import EventCollection
import math
from scipy import integrate
from scipy import constants

K_b = 1.380649e-23  # Boltzmanns Constant
a_n = 6.022e23  # Avogrados Number
T = 10  # Kelvin
mass = 28
a = 10e-5

def mantleEvolution(tMax, dt) :
    N = int(tMax/dt) + 1
    timeArray = np.zeros(N)
    numberArray = np.zeros(N)
    t = 0
    n_m = 1 #per cubic centimeter
    n_s = 0
    i = 0
    while i < N:
        m_m = float(mass) / a_n
        V_m = math.sqrt((8.0 * K_b * float(T)) / (math.pi * float(m_m)))
    R = math.pi * pow(float(a), 2) * n_m * V_m
    ρ_d = 3  # g/cm^3
    μ_g = 2 * 1.67e-23
    m_d = (4 / 3) * math.pi * pow(float(a), 3)
    n_g = 10e4
    d_g = 0.01  # 1% mass density of gas
    n_d = (d_g * n_g * μ_g) / (m_d * float(ρ_d))
    n_s = n_s + R * dt
    n_m = n_m - R * n_d * dt
    timeArray[i] = t
    numberArray[i] = n_m
    t = t + dt
    i = i + 1
return [timeArray, numberArray, n_s]

timeArray, numberArray, n_s = mantleEvolution(5.0 * 10**15, 5.0 * 10**9)

fig = plt.figure(figsize=(6.4, 6.4))
ax1 = plt.subplot(111)
ax1.plot(timeArray, numberArray)
ax1.set_xlabel('time, seconds', fontsize=20)
ax1.set_ylabel('number', fontsize=20)
plt.setp(ax1.get_xticklabels(), fontsize=16)
plt.setp(ax1.get_yticklabels(), fontsize=16)
fig.subplots_adjust(left=.18)
plt.savefig('mantleEvolution.pdf')

原运行结果 #

需求 #

绘制第二条曲线，使其增长速率与现有曲线的下降速率完全一致，效果参考下图：

修改后的代码 #

import scipy
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.collections import EventCollection
import math
from scipy import integrate
from scipy import constants

K_b = 1.380649e-23  # 玻尔兹曼常数
a_n = 6.022e23  # 阿伏伽德罗常数
T = 10  # 开尔文
mass = 28
a = 10e-5

def mantleEvolution(tMax, dt) :
    N = int(tMax/dt) + 1
    timeArray = np.zeros(N)
    numberArray = np.zeros(N)
    # 新增数组存储每个时间点的n_s值
    numberArray_s = np.zeros(N)
    t = 0
    n_m = 1  # 每立方厘米
    n_s = 0
    i = 0
    while i < N:
        m_m = float(mass) / a_n
        V_m = math.sqrt((8.0 * K_b * float(T)) / (math.pi * float(m_m)))
        R = math.pi * pow(float(a), 2) * n_m * V_m
        ρ_d = 3  # g/cm^3
        μ_g = 2 * 1.67e-23
        m_d = (4 / 3) * math.pi * pow(float(a), 3)
        n_g = 10e4
        d_g = 0.01  # 气体质量密度的1%
        n_d = (d_g * n_g * μ_g) / (m_d * float(ρ_d))
        # 让n_s的增量等于n_m的减量，实现速率完全反向匹配
        n_s = n_s + R * n_d * dt
        n_m = n_m - R * n_d * dt
        timeArray[i] = t
        numberArray[i] = n_m
        numberArray_s[i] = n_s
        t = t + dt
        i = i + 1
    return [timeArray, numberArray, numberArray_s]

timeArray, numberArray, numberArray_s = mantleEvolution(5.0 * 10**15, 5.0 * 10**9)

fig = plt.figure(figsize=(6.4, 6.4))
ax1 = plt.subplot(111)
ax1.plot(timeArray, numberArray, label='地幔粒子数 (n_m)')
# 绘制反向增长的第二条曲线
ax1.plot(timeArray, numberArray_s, label='增长粒子数 (n_s)')
ax1.set_xlabel('时间（秒）', fontsize=20)
ax1.set_ylabel('数量', fontsize=20)
plt.setp(ax1.get_xticklabels(), fontsize=16)
plt.setp(ax1.get_yticklabels(), fontsize=16)
# 添加图例区分曲线
ax1.legend(fontsize=16)
fig.subplots_adjust(left=.18)
plt.savefig('mantleEvolution.pdf')
plt.show()

关键修改点 #

• 修正增长逻辑：将n_s的更新公式从n_s += R*dt改为n_s += R*n_d*dt，确保其增长速率与n_m的下降速率绝对值完全相等，实现反向同步变化。
• 记录时间序列：新增numberArray_s数组，保存每个时间点的n_s数值，用于绘制随时间变化的曲线。
• 优化可视化：添加图例、中文标签，提升图表可读性。

内容的提问来源于stack exchange，提问作者Santiago