``````

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 [ ]:

``````