In [1]:
import sympy as smp
x, y, t = smp.symbols('x, y, t')
import numpy as np
pi = np.pi
Define an abribrary solution.
In [2]:
#u = smp.sin(0.5*pi*t) * smp.cos(3*pi*(x-1.0)) * smp.cos(7*pi*y)
u = smp.cos(0.5*pi*t) * smp.cos(3*pi*(x-1.0)) * smp.cos(7*pi*y)
#u = (1-smp.cos(0.5*pi*t)) * smp.cos(3*pi*(x-1.0)) * smp.cos(7*pi*y)
#u = smp.log(1+t) * smp.cos(3*pi*(x-1.0)) * smp.cos(7*pi*y)
Define matching right-hand side for heat equation $u_t - \Delta u = f$.
In [3]:
f = smp.diff(u,t) - smp.diff(smp.diff(u, x), x) - smp.diff(smp.diff(u, y), y)
f = smp.expand(f)
Define forward Euler scheme:
In [5]:
def forward_euler(u0, dt, f0):
u1 = u0 + dt * (smp.diff(u0, x, 2) + smp.diff(u0, y, 2)) + dt * f0
#u1 = u0
return u1
Perform one step with the above solution.
In [6]:
Dt = [0.5**k for k in range(10)]
Err = []
for dt in Dt:
u0 = u.subs({t:0.0})
f0 = f.subs({t:0.0})
u_approx = forward_euler(u0, dt, f0)
u1 = u.subs({t:dt})
err = u1 - u_approx
# Compute L^2 norm over unit square [0,1]^2.
err_norm = smp.sqrt(smp.integrate(err**2, (x, 0.0, 1.0), (y, 0.0, 1.0)))
Err.append(err_norm)
Plot error norm.
In [9]:
from matplotlib import pyplot as pp
# Show orders 1,2,3 for comparison
for o in [1,2,3,4]:
pp.loglog([Dt[0], Dt[-1]], [Err[0], Err[0]*(Dt[-1]/Dt[0])**o], color='0.5')
pp.loglog(Dt, Err,'o-')
pp.show()
In [66]:
In [ ]: