使用python integrator检查收敛性

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