I 'm trying to solve the 4 types of heat equations homogeneous linear, non-homogeneous and non linear. This is the code that I have written for linear homogeneous case for beta = -1 and N of your choice. I have used hat functions with x step size = H and time I 'm wishing to take $H^2$. The problem is that if I 'm using scipy's runge kutta's method with automatic step size it's giving me this graph enter image description here which suddenly gives high error at 256 basis but for basis less than this I 'm getting the error = c $H^p$ log taken relationship. So I can't say it's converging but it's taking too much time for further computation. If I run the algo_2 in the code which I have copied from geeksforgeeks or even if I write own as given on wikipedia I get infinite in just 4 time loop. Can someone please tell me what's wrong here? Equation is $u_t + \beta u_{xx} = f(x, t)$ and the weak derivative form uses $\phi_j$ as the hat functions in [0, 1]. I 'm trying to compute the values at t=1. I am also attaching the link to github. https://github.com/XtremeGood/Maths-SOP/tree/corrected_boundary_conditions/working%201 Please check the corrected_boundary_conditions branch's working1 folder for the graph image and files if it doesn't open.
class Heat:
#A *alpha - beta(D *) +beta(E(1 - 0)alpha) = 0
def __init__(self, N=10):
self.N = N
self.H = 1/N
self.domain = np.linspace(0, 1, N+1)
self.beta = -1
self.t = 1
#Set initial conditions here
self.alpha = np.sin(np.pi* self.domain)
#Set solution here
self.exact = lambda t:math.exp(t) * self.alpha
#Define boundary conditions
self.b_c_start = lambda t:0
self.b_c_end = lambda t:0
#Define f
def f(k):
if k == 0:
return (-(math.pi*self.H - math.sin(math.pi*self.H))*(math.pi**2*self.beta - 1)/(math.pi**2*self.H))
elif k == N:
return ((-math.pi*self.H + math.sin(math.pi*self.H))*(math.pi**2*self.beta - 1)/(math.pi**2*self.H))
else:
return (2*(math.pi**2*self.beta - 1)*(math.cos(math.pi*self.H) - 1)*math.sin(math.pi*self.H*k)/(math.pi**2*self.H))
def c_dash_transform(l):
l[0] /= self.A[0]
for i in range(1, self.N):
l[i] = (l[i])/ (self.A[i]-(self.H/6)* l[i-1])
return l
self.A = self.H* np.array([1/3] + [2/3]*(N-1) + [1/3])
# self.A = (self.H) * sp.csc_matrix(diags([ [1/6 for i in range(N)], \
# [1/3] + [2/3]*(N-1) + [1/3], \
# [1/6 for i in range(N)] ], \
# [1, 0, -1]), dtype=np.float32)
self.c_dash = c_dash_transform(self.H *np.array( [1/6]*(N)))
#-beta*
self.D = N * sp.csr_matrix(diags([ [-1 for i in range(N)], \
[1] + [2]*(N-1) + [1], \
[-1 for i in range(N)]], \
[1, 0, -1] ), dtype=np.float32)
self.D[0, 0] -= N
self.D[0, 1] -= -N
self.D[N, N-1] -= -N
self.D[N, N] -= N
self.F = np.array([f(k) for k in prange(N+1)]).reshape(-1, 1) #incomplete without exp(t)
def error(self):
return np.linalg.norm(self.exact(self.t) - self.approx, 2)
def func (self, t, y):
bc_s = [self.b_c_start(t),self.b_c_end(t)]
def apply_bc(y):
y[0, 0], y[-1, 0] = bc_s
return y
def thomas_algorithm(d ):
d[0] /= self.A[0]
for i in range(1, self.N + 1):
d[i, 0] = (d[i, 0] - (self.H/6)*d[i-1])/(self.A[i]-(self.H/6)*self.c_dash[i-1])
for i in range(self.N-1, -1, -1): # Last entry is left unchanged ;see that i+1
d[i, 0] -= self.c_dash[i]*d[i+1]
return d
y = apply_bc(y)
d = self.beta*self.D*y + math.exp(t)*self.F # A alpha' = d
d = thomas_algorithm(d) # alpha' = d
d = apply_bc(d)
return d
def algo_1(self):
self.algo = solve_ivp(self.func,t_span=(0, self.t + .1), y0=self.alpha, \
t_eval= np.array([self.t]), vectorized=True, method="LSODA")
self.algo.y[0, 0], self.algo.y[-1, 0] =[self.b_c_start(self.t),self.b_c_end(self.t)]
self.approx = self.algo.y.reshape(-1, )
def algo_2(self):
# Finds value of y for a given x using step size h
# and initial value y0 at x0.
def rungeKutta(y):
t = 0
while t != self.t:
"Apply Runge Kutta Formulas to find next value of y"
k1 = self.H**2 * self.func(t, y)
k2 = self.H**2 * self.func(t + 0.5 * self.H**2, y + 0.5 * k1)
k3 = self.H**2 * self.func(t + 0.5 * self.H**2, y + 0.5 * k2)
k4 = self.H**2 * self.func(t + self.H**2, y + k3)
# Update next value of y
y = y + (1.0 / 6.0)*(k1 + 2 * k2 + 2 * k3 + k4)
# Update next value of x
t += self.H**2
return y
self.algo = rungeKutta(self.alpha.reshape(-1, 1))
self.algo[0, 0], self.algo[-1, 0] = [self.b_c_start(self.t),self.b_c_end(self.t)]
self.approx = self.algo.reshape(-1, )