8.5 求微分方程组(竖直加热板的自然对流)的数值解。(d^3f)/(dm^3)+3f(d^f)/(dm^2)-2((df)/(dm))^2+T=0,(d^2T)/(dm^2)+2.1f(dT)/(dm)=0,已知当m=0时,f=0,(df)/(dm)=0,(d^2f)/(dm^2)=0.68,T=1,(dT)/(dm)=-0.5。要求在[0,10]区间上,画出f(a)、T(a)的解曲线

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp

def model(t, y):
f, df_dm, d2f_dm2, T, dT_dm = y
d3f_dm3 = -3fd2f_dm2 + 2(df_dm)**2 - T
d2T_dm2 = -2.1
f*dT_dm
return [df_dm, d2f_dm2, d3f_dm3, dT_dm, d2T_dm2]

y0 = [0, 0, 0.68, 1, -0.5]

t_span = (0, 10)
t_eval = np.linspace(t_span[0], t_span[1], 1000)

sol = solve_ivp(model, t_span, y0, t_eval=t_eval, method='RK45')

f = sol.y[0]
T = sol.y[3]

plt.figure(figsize=(12, 6))

plt.subplot(2, 1, 1)
plt.plot(sol.t, f, label='f(m)')
plt.xlabel('m')
plt.ylabel('f(m)')
plt.title('Solution for f(m)')
plt.grid(True)

plt.subplot(2, 1, 2)
plt.plot(sol.t, T, label='T(m)', color='orange')
plt.xlabel('m')
plt.ylabel('T(m)')
plt.title('Solution for T(m)')
plt.grid(True)

plt.tight_layout()
plt.show()

print("学号后四位:04")

posted on   黄元元  阅读(26)  评论(0编辑  收藏  举报

相关博文:
阅读排行:
· winform 绘制太阳,地球,月球 运作规律
· AI与.NET技术实操系列(五):向量存储与相似性搜索在 .NET 中的实现
· 超详细:普通电脑也行Windows部署deepseek R1训练数据并当服务器共享给他人
· 【硬核科普】Trae如何「偷看」你的代码?零基础破解AI编程运行原理
· 上周热点回顾(3.3-3.9)
< 2025年3月 >
23 24 25 26 27 28 1
2 3 4 5 6 7 8
9 10 11 12 13 14 15
16 17 18 19 20 21 22
23 24 25 26 27 28 29
30 31 1 2 3 4 5

统计

点击右上角即可分享
微信分享提示