I am looking to integrate the difference between my numerical and exact solution to the heat equation though I am not sure what would be the best to way to tackle this. Is there a specific integrator that would allow me to do this ? I hope to integrate it wrt to $x$.
I have the following code so far:
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)
What I am looking to integrate is the variable diff (the squared difference). Any suggestions are welcomed, thanks.
Edit: I would like to use RK45 method however I am not sure what should my y0 be, I have tried x1 = integrate.RK45(diff,diff[0],0,t_end) but get the following output error:
raise ValueError("`y0` must be 1-dimensional.")
ValueError: `y0` must be 1-dimensional.
