一、题目
二、代码
from scipy.integrate import quad
import numpy as np
# 函数f(u)=1/2*u**2,故f`(u)=u
def f_positive(upp_value):
# 积分f+(u)中需要使用的函数
is_over_zero = int(upp_value > 0)
return is_over_zero * upp_value
def f_negative(upp_value):
# 积分f-(u)中需要使用的函数
is_over_zero = int(upp_value > 0)
return (1 - is_over_zero) * upp_value
def cal_fu(upp_value, form='+'):
"""
:param upp_value:积分上限
:param form: 可选值列表[‘+’,‘-’],决定返回f+(u)还是f-(u)
:return: quad_r: f+(u)、f-(u)的值
"""
quad_r = 0
if form == "+":
quad_r = quad(f_positive, 0, upp_value)[0]
elif form == '-':
quad_r = quad(f_negative, 0, upp_value)[0]
return quad_r
def cal_next_step(pt, px):
"""
该空间点在下一时间层的值
:param pt: 当前点在网格的空间位置
:param px: 当前点在网格的时间位置
:param grid: 网格点上的值
:return:
"""
u_j_n = grid_value[pt][px] # 计算U(j,n)
u_jplus1_n = grid_value[pt][px+1] # 计算U(j+1,n)
u_jminus1_n = grid_value[pt][px-1] # 计算U(j-1,n)
temp_minus = cal_fu(u_jplus1_n, form='-') - cal_fu(u_j_n, form='-')
temp_plus = cal_fu(u_j_n) + cal_fu(u_jminus1_n)
u_j_nplus1 = u_j_n - grid_ratio * (temp_minus + temp_plus)
grid_value[pt+1][px] = u_j_nplus1
if __name__ == "__main__":
x_range = [-2, 2] # 空间范围
t_range = [0, 0.9] # 时间范围
delta_x = 0.1 # 空间步长
delta_t = 0.01 # 时间步长
grid_ratio = delta_t / delta_x # 网格比
grid_x = int((x_range[1] - x_range[0]) / delta_x) + 1 # 空间网格点数,此例中为41
grid_t = int((t_range[1] - t_range[0]) / delta_t) + 1 # 时间网格点数,此例中为91
# 考虑用列表grid_value来存储Ujn[[t=0.01],...,[t=0.9]]
grid_value = np.zeros((grid_t, grid_x)) # 行代表某个时间、列代表某个空间
# 将初始值t=0添加到grid_value中,即初始条件
for i in range(grid_x):
x_current = x_range[0] + delta_x * i
if x_current > 0:
grid_value[0][i] = 1
else:
grid_value[0][i] = -1
# 将每一个时间层上的左右边界赋固定值
grid_value[:, grid_x-1] = 1 # 右边界为1
grid_value[:, 0] = -1 # 左边界为-1
# 开始计算,时间上索引从0算到89,空间上索引从1算到39
# 假设右边界必定收敛
for i in range(0, grid_t-1):
for j in range(1, grid_x-1):
cal_next_step(i, j)
grid_value[i, -1] = grid_value[i, -2]
grid_value[-1, -1] = grid_value[-1, -2]
# 仅将最后一个时间层的网格点数据保存到"2.txt"中
np.savetxt('2.txt', grid_value[-1, :], fmt='%0.8f')
三、运行结果
绘图如下: