使用 python 积分器检查收敛

Question

我希望整合热方程的数值解和精确解之间的差异，尽管我不确定解决这个问题的最佳方法是什么。是否有特定的集成商可以让我这样做？ 我希望将它集成到 $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.

Answer 1

通过集成，您的意思是要找到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