Heat equation inversion
The heat equation is another classic example of an inverse problem. The following example shows how you can estimate time dependent temperature distribution on a 1d domain using several point measurements when heating is localized to a few fixed points. This is done with the help of sparse matrices and regularized least squares. Regularization of the second derivative used used to ensure smoothness of the reconstructed heating.
A generalization of this to 2d and 3d (and with advection added to the model) has applications for example in oceanography, where the inverse problem is to estimate the methane release from localized sources in the sea floor.
#!/usr/bin/env python
#
# Inverse problems with the 1d heat equation
# (c) 2017 Juha Vierinen
#
import numpy as n
import matplotlib.pyplot as plt
from scipy.sparse import coo_matrix
import scipy.sparse.linalg as s
import lsqlin
def create_theory_matrix(n_t=200,n_x=200,
dt=2.0,dx=1.0,C=0.5,
alpha=0,
alpha_q=0,
alpha_u=0,
q_true=None,q_x=100, # measured heating at location x
u_meas=None,u_x=100): # meuasred u at location x
#
# n_t = number of points in time
# n_x = number of points in space
# dx = spatial grid spaceing
# dt = time spacing
# alpha = L_2 regularization
# parameters for u and
# n_t parameters for heat source q at known location (center of domain)
#
print(q_x)
n_par=n_t*n_x + n_t*len(q_x)
qi0=n_t*n_x
# diffusion coefficient (assumed to be known)
C = n.zeros([n_t,n_x])
C[:,:]=0.5
# number of unknowns
n_col=n_par
row_idx=0
d=[]
rows=[]
cols=[]
data=[]
row_idx=0
q_x=n.array(q_x)
# heat equation part du/dx - C d^2u/dx^2 - q = 0
alpha_theory=1.0
for i in range(n_t-1):
for j in range(1,n_x-1):
# du/dt
rows.append(row_idx)
cols.append(i*n_x+j)
data.append(alpha_theory*-1.0/dt)
rows.append(row_idx)
cols.append((i+1)*n_x+j)
data.append(alpha_theory*1.0/dt)
# -C*d^2u/dx^2
rows.append(row_idx)
cols.append(i*n_x+j-1)
data.append(alpha_theory*-C[i,j]*1.0/(dx**2.0))
rows.append(row_idx)
cols.append(i*n_x+j)
data.append(alpha_theory*C[i,j]*2.0/(dx**2.0))
rows.append(row_idx)
cols.append(i*n_x+j+1)
data.append(alpha_theory*-C[i,j]*1.0/(dx**2.0))
# specific location at center where heating occurs
if j in q_x:
q_idx=n.where(j==q_x)[0][0]
rows.append(row_idx)
cols.append(qi0+q_idx*n_t+i)
data.append(alpha_theory*-1.0)
d.append(0.0)
row_idx+=1
# initial boundary condition of 0
ut=n.repeat(0,n_x)
ux=n.arange(n_x)
uv=n.zeros(n_x)
for mi in range(len(ut)):
i=ut[mi]
j=ux[mi]
rows.append(row_idx)
cols.append(i*n_x+j)
data.append(1.0)
d.append(uv[mi])
row_idx+=1
# regularization of heat field
if alpha > 0:
for i in range(n_t-2):
for j in range(n_x):
rows.append(row_idx)
cols.append(i*n_x+j)
data.append(1.0*alpha)
rows.append(row_idx)
cols.append((i+1)*n_x+j)
data.append(-2.0*alpha)
rows.append(row_idx)
cols.append((i+2)*n_x+j)
data.append(1.0*alpha)
d.append(0.0)
row_idx+=1
# regularization of heat sources
if alpha_q > 0:
for j in range(len(q_x)):
for i in range(n_t-2):
rows.append(row_idx)
cols.append(qi0+j*n_t+i)
data.append(1.0*alpha_q)
rows.append(row_idx)
cols.append(qi0+j*n_t+i+1)
data.append(-2.0*alpha_q)
rows.append(row_idx)
cols.append(qi0+j*n_t+i+2)
data.append(1.0*alpha)
d.append(0.0)
row_idx+=1
# heating boundar conditions q=0 at t=0 and q=0 at t=n_t
for j in range(len(q_x)):
rows.append(row_idx)
cols.append(qi0+j*n_t)
data.append(1.0)
d.append(0.0)
row_idx+=1
rows.append(row_idx)
cols.append(qi0+j*n_t+n_t-1)
data.append(1.0)
d.append(0.0)
row_idx+=1
# heating measurements (used for forward model)
if q_true != None:
for j in range(len(q_x)):
for i in range(n_t):
rows.append(row_idx)
cols.append(qi0+j*n_t+i)
data.append(1.0)
d.append(q_true[j][i])
row_idx+=1
if u_meas != None:
for j in range(len(u_x)):
for i in range(n_t):
rows.append(row_idx)
cols.append(i*n_x+u_x[j]) # at time i, location u_x
data.append(alpha_u)
d.append(u_meas[j][i]*alpha_u)
row_idx+=1
d=n.array(d)
n_row=row_idx
n_col=n_par
# create a sparse theory matrix
G=coo_matrix((data, (rows, cols)), shape=(n_row, n_col), dtype=n.float32)
return(G,d,qi0)
n_t=100 # number of time steps
n_x=200 # number of grid points
# true heating (gaussian shape)
q_true0=n.exp(-0.5*(n.linspace(-5,5,num=n_t)+1.0)**2.0/0.1)
q_true0=q_true0/n.max(q_true0)
q_true1=n.exp(-0.5*(n.linspace(-5,5,num=n_t)-1.0)**2.0/0.1)
q_true1=q_true1/n.max(q_true1)
q_true=[q_true0,q_true1]
q_x=[82,140] # heat source locations
u_x=[50,90,150] # temperature measurement location
C=0.2 # diffusion coefficient
# solve heat equation, no regularization needed, because the problem is fully determined
# as we provide heating q.
G,d,qi0=create_theory_matrix(n_t=n_t,n_x=n_x,q_true=q_true,q_x=q_x,alpha=1e-6,alpha_q=0,C=C)
print(G.shape)
print(d.shape)
m0=s.lsqr(G,d)[0]
m0u=m0[0:qi0]
q_hat0=m0[qi0+n.arange(n_t)]
q_hat1=m0[qi0+n_t+n.arange(n_t)]
m0u.shape=(n_t,n_x)
# simulate measurement, add noise
u_meas_0=m0u[:,u_x[0]]+n.random.randn(n_t)*0.02
u_meas_1=m0u[:,u_x[1]]+n.random.randn(n_t)*0.02
u_meas_2=m0u[:,u_x[2]]+n.random.randn(n_t)*0.02
u_meas=[u_meas_0,u_meas_1,u_meas_2]
# now create a theory matrix with q unknown
# in order to estimate q from a measurement of u from only one location
#
# 1.0/alpha_u is "measurement error standard deviation"
# the larger alpha_u, the more weight we put on measurements and
# less weight on the solution satisfying the heat equation
q_x=[42,82,140] # heat source locations
G,d,qi0=create_theory_matrix(n_t=n_t,
n_x=n_x,
u_meas=u_meas,
q_x=q_x,
u_x=u_x,
alpha=0.01,alpha_q=0.01,
alpha_u=0.1,C=C)
print(G.shape)
print(d.shape)
# sparse iterative least squares
# it would be advantageous to add a positivity constraint to u and q, but maybe later
m0=s.lsqr(G,d)[0]
m0u_hat=m0[0:qi0]
m0u_hat.shape=(n_t,n_x)
q_hat0=m0[qi0+n.arange(n_t)]
q_hat1=m0[qi0+n_t+n.arange(n_t)]
q_hat2=m0[qi0+2*n_t+n.arange(n_t)]
#
# Plot stuff
#
plt.subplot(231)
plt.plot(q_true[0])
plt.plot(q_true[1])
plt.title("True heating")
plt.xlabel("Time")
plt.ylim([-0.5,1.5])
#plt.title("$q_{\mathrm{true}}$ at x",q_x)
plt.subplot(232)
plt.pcolormesh(n.transpose(m0u))
plt.ylabel("x")
plt.xlabel("Time")
plt.axhline(q_x[1],color="red")
plt.axhline(q_x[2],color="red")
plt.axhline(u_x[0],color="white")
plt.axhline(u_x[1],color="white")
plt.axhline(u_x[2],color="white")
plt.title("Temperature (true)")
plt.colorbar()
plt.subplot(233)
plt.plot(u_meas_0,label="Measurement")
plt.plot(u_meas_1,label="Measurement")
plt.plot(u_meas_2,label="Measurement")
plt.title("Measurements")
plt.legend(loc=2)
plt.xlabel("Time")
plt.ylabel("Temperature")
#plt.title("Measured temperature at x ",u_x)
plt.subplot(234)
plt.plot(q_hat0,label="Estimate",color="black")
plt.plot(q_hat1,label="Estimate",color="black")
plt.plot(q_hat2,label="Estimate",color="red")
plt.plot(q_true[0],label="True")
plt.plot(q_true[1],label="True")
plt.legend()
plt.ylim([-0.5,1.5])
plt.title("Estimated heating")
#plt.title("$\hat{q}$ at x",q_x)
plt.subplot(235)
plt.pcolormesh(n.transpose(m0u_hat),vmin=n.min(m0u),vmax=n.max(m0u))
plt.axhline(u_x[0],color="white")
plt.axhline(u_x[1],color="white")
plt.axhline(u_x[2],color="white")
plt.axhline(q_x[0],color="red")
plt.axhline(q_x[1],color="red")
plt.axhline(q_x[2],color="red")
plt.ylabel("x")
plt.xlabel("Time")
plt.title("Estimated temperature")
plt.colorbar()
plt.subplot(236)
plt.plot(m0u_hat[:,u_x[0]],label="Estimated")
plt.plot(m0u_hat[:,u_x[1]],label="Estimated")
plt.plot(m0u_hat[:,u_x[2]],label="Estimated")
plt.plot(m0u[:,u_x[0]],label="True")
plt.plot(m0u[:,u_x[1]],label="True")
plt.plot(m0u[:,u_x[2]],label="True")
plt.legend(loc=2)
#plt.title("Estimated temperature at x ",u_x)
plt.xlabel("Time")
plt.show()
Comments
Post a Comment