python - 使用 python 积分器检查收敛
问题描述
我希望将热方程的数值解和精确解之间的差异整合起来,尽管我不确定解决这个问题的最佳方法是什么。是否有特定的集成商可以让我这样做?我希望将它集成到 $x$。
到目前为止,我有以下代码:
import numpy as np
from matplotlib import pyplot
from mpl_toolkits.mplot3d import Axes3D
from scipy import linalg
import matplotlib.pyplot as plt
import math
def initial_data(x):
y = np.zeros_like(x)
for i in range(len(x)):
if (x[i] < 0.25):
y[i] = 0.0
elif (x[i] < 0.5):
y[i] = 4.0 * (x[i] - 0.25)
elif (x[i] < 0.75):
y[i] = 4.0 * (0.75 - x[i])
else:
y[i] = 0.0
return y
def heat_exact(x, t, kmax = 150):
"""Exact solution from separation of variables"""
yexact = np.zeros_like(x)
for k in range(1, kmax):
d = -8.0*
(np.sin(k*np.pi/4.0)-2.0*np.sin(k*np.pi/2.0)+np.sin(3.0*k*np.pi/4.0))/((np.pi*k)**2)
yexact += d*np.exp(-(k*np.pi)**2*t)*np.sin(k*np.pi*x)
return yexact
def ftcs_heat(y, ynew, s):
ynew[1:-1] = (1 - 2.0 * s) * y[1:-1] + s * (y[2:] + y[:-2])
# Trivial boundary conditions
ynew[0] = 0.0
ynew[-1] = 0.0
Nx = 198
h = 1.0 / (Nx + 1.0)
t_end = 0.25
s = 1.0 / 3.0 # s = delta / h^2
delta = s * h**2
Nt = int(t_end / delta)+1
x = np.linspace(0.0, 1.0, Nx+2)
y = initial_data(x)
ynew = np.zeros_like(y)
for n in range(Nt):
ftcs_heat(y, ynew, s)
y = ynew
fig = plt.figure(figsize=(8,6))
ax = fig.add_subplot(111)
x_exact = np.linspace(0.0, 1.0, 200)
ax.plot(x, y, 'kx', label = 'FTCS')
ax.plot(x_exact, heat_exact(x_exact, t_end), 'b-', label='Exact solution')
ax.legend()
plt.show()
diff = (y - heat_exact(x_exact,t_end))**2 # squared difference between numerical and exact solutions
x1 = np.trapz(diff, x=x) #(works!)
import scipy.integrate as integrate
x1 = integrate.RK45(diff,diff[0],0,t_end) #(preferred but does not work)
我要整合的是变量diff
(平方差)。欢迎任何建议,谢谢。
编辑:我想使用 RK45 方法但是我不确定我的 y0 应该是什么,我已经尝试过x1 = integrate.RK45(diff,diff[0],0,t_end)
但得到以下输出错误:
raise ValueError("`y0` must be 1-dimensional.")
ValueError: `y0` must be 1-dimensional.
解决方案
通过集成,您的意思是您想要找到 和 之间的y
区域heat_exact
?或者您想知道它们在特定限制内是否相同?后者可以通过 找到numpy.isclose
。前者你可以使用 numpy 内置的几个集成函数。
例如:
np.trapz(diff, x=x)
哦,最后一行不应该是diff = (y - heat_exact(x_exact,t_end))**2
吗?我对此的整合diff
给出了 8.32E-12,从你给我的情节来看,这看起来是正确的。
还可以查看scipy.integrate
推荐阅读
- javascript - 如何在下一个输入元素中获取选择选项数据
- mysql - MySQL - 比较同一张表中的不同键和值
- firebase - 将 Chrome 扩展程序连接到 Firebase RTB 时出错:“firebase.database 不是函数”
- android - JADX 从混淆的应用程序中获取类名
- android - Android 构建停止使用 ClassNotFoundException
- typo3-9.x - 如何始终从 TYPO3 后端的图像中激活 alt 标签?
- python - 使用 Python 将多个 xml 文件/链接转换为 JSON?
- java - 需要从另一个VM访问API(https://localhost:8080/er/heartbeats)
- python - 如何编译单个python脚本(不是exe)?
- python - 使用 MongoDB PyMongo 进行计数的聚合函数