# Is this the correct way for solving coupled 1d PDEs using finite difference methods?

I am trying to solve the following coupled PDEs:

$$C_e\frac{\partial u(x,t)}{\partial t} = k_{ed}\frac{\partial^2u(x,t)}{\partial^2x} - G_{el}(u(x,t) - v(x,t)) + S(x,t)$$ $$C_l\frac{\partial v(x,t)}{\partial t} = k_{ed}\frac{\partial^2v(x,t)}{\partial^2x} + G_{el}(u(x,t) - v(x,t))$$

where $$k_{ed}$$, $$G_{el}$$ are constants and $$C_e = \gamma u(x,t)$$ and $$C_l$$ are specific heats and $$S(x,t) = I_0 exp(-\frac{x}{\delta}-\frac{(t - 2t_p)^2}{t_p^2})$$ is the laser pulse source term.

I use finite difference methods to solve the above equations as follows:

$$u^{f+1}_i = u^f_i + \frac{k_{ed}dt}{\Delta x^2}(u^f_{i+1} -2u^f_i + u^f_{i-1}) + dt(-G_{el}(u^f_i - v^f_i) + S^f_i)$$ and $$v^{f+1}_i = v^f_i + \frac{k_{ed}dt}{\Delta x^2}(v^f_{i+1} -2v^f_i + v^f_{i-1}) + dt(G_{el}(u^f_i - v^f_i))$$

Where $$(f,i)$$ are mesh in time and space as $$F\Delta t = T_{total}, f = 0,1,2 ... F$$ and $$N\Delta x = L, i = 0,1,2,... N$$ where $$T_{total}$$ and $$L$$ are the total integration time and length respectively. I have implemented the above in python:

import numpy as np
import matplotlib.pyplot as plt
# sample parameter §§ Gold
la0 = 429 # conductivity in W/mK
gma = 62.8 # thermal constant J/m^3K^2
Cl = 2.6*10**6 # phonon specific heat  in J/m^3K
Gel = 3.5*10**16 # lelectron phonon coupling constant
tau_e = 0.04 # e relaxation t const
tau_l = 0.6 # l relaxation t const
Tm = 300
L = 100 # sample thickness in nm
# Laser pulse parameters
T = 6

R = 0.93 # reflectivity
I0 = 100 # fulence in J/m^2
tp = 0.1 # in ps
zz0 = 15.3 # penetration depth in nm

#F = 0.5

def las(x,t):
pt = 0.93*((1-R)/(zz0*tp))*I0*np.exp(-(x/zz0)-2.772*((t-2*tp)**2/tp**2))
return pt

def I(x):
tin = Tm
return tin

a = la0 # Ked

def solver_FE(I, a, las, L, T):
"""

"""
import time;  t0 = time.clock()  # For measuring the CPU time

dt = 0.001
Nt = int(round(T/float(dt)))
t = np.linspace(0, Nt*dt, Nt+1)   # Mesh points in time

dx = 1
Nx = int(round(L/dx))
x = np.linspace(0, L, Nx+1)       # Mesh points in space

F = a*(dt/dx**2)
u   = np.zeros(Nx+1)
u_n = np.zeros(Nx+1)
v   = np.zeros(Nx+1)
v_n = np.zeros(Nx+1)

# Set initial condition u(x,0) = I(x)
for i in range(0, Nx+1):
u_n[i] = I(x[i])
v_n[i] = I(x[i])

for n in range(0, Nt):
# Compute u, v  at inner mesh points
for i in range(1, Nx):
u[i] = u_n[i] + (-F*(u_n[i-1] - 2*u_n[i] + u_n[i+1]) + dt*(-Gel*(u_n[i] - v_n[i]) + las(x[i], t[n])))/(gma*u_n[i])
v[i] = v_n[i] + (-F*(v_n[i-1] - 2*v_n[i] + v_n[i+1]) + dt*(Gel*(u_n[i] - v_n[i])))/(Cl)

# Insert boundary conditions
u = 0;  u[Nx] = 0
v = 0;  v[Nx] = 0

# this step is to save each value of v and u in each time step before updating it, such that I can plot u and v with respect to time.
ut = np.array([])
vt = np.array([])
for m in range(0, Nt):
ut = np.append(ut,u)
vt = np.append(vt,v)

# Switch variables before next step

u_n, u = u, u_n
v_n, v = v, v_n

t1 = time.clock()
return u_n, v_n, x, t, t1-t0  # u_n holds latest u

u, v, x, t, cpu = solver_FE(I, a, las, L, T)

#fig = plt.figure(1)
##t = np.linspace(0, T, h)
#plt.plot(t,u,'r',label=r'$$T_e$$')
#plt.plot(t,v,'b',label=r'$$T_l$$')
##plt.plot(t,sol0[:,2],'g',label=r'$$T_s$$')
#plt.ylabel('$$Temperature$$',fontsize=20)
#plt.xlabel('Delay',fontsize=20)
#plt.legend(loc='best')
#plt.show()
##
fig = plt.figure(2)
#t = np.linspace(0, T, h)
plt.plot(x,u,'r',label=r'$$T_e$$')
plt.plot(x,v,'g',label=r'$$T_l$$')
#plt.plot(t,sol0[:,2],'g',label=r'$$T_s$$')
plt.ylabel('$$Temperature$$',fontsize=20)
plt.xlabel('length',fontsize=20)
plt.legend(loc='best')
plt.show()

fig = plt.figure(3)
#t = np.linspace(0, T, h)
plt.plot(t,las(0,t),'r',label=r'$$T_e$$')
#plt.plot(x,v,'g',label=r'$$T_l$$')
#plt.plot(t,sol0[:,2],'g',label=r'$$T_s$$')
plt.ylabel('$$Intensity$$',fontsize=20)
plt.xlabel('Delay',fontsize=20)
plt.legend(loc='best')
plt.show()


while running it I get the following errors:

RuntimeWarning: overflow encountered in double_scalars
u[i] = u_n[i] + F*(u_n[i-1] - 2*u_n[i] + u_n[i+1]) + dt*(-Gel*(u_n[i] - v_n[i]) + las(x[i], t[n]))
C:/Users/jayas/OneDrive/Documents/Python Scripts/llg/diff_ttm_v_test.py:90: RuntimeWarning: overflow encountered in double_scalars
v[i] = v_n[i] + F*(v_n[i-1] - 2*v_n[i] + v_n[i+1]) + dt*(Gel*(u_n[i] - v_n[i]))
C:/Users/jayas/OneDrive/Documents/Python Scripts/llg/diff_ttm_v_test.py:89: RuntimeWarning: invalid value encountered in double_scalars
u[i] = u_n[i] + F*(u_n[i-1] - 2*u_n[i] + u_n[i+1]) + dt*(-Gel*(u_n[i] - v_n[i]) + las(x[i], t[n]))
C:/Users/jayas/OneDrive/Documents/Python Scripts/llg/diff_ttm_v_test.py:90: RuntimeWarning: invalid value encountered in double_scalars
v[i] = v_n[i] + F*(v_n[i-1] - 2*v_n[i] + v_n[i+1]) + dt*(Gel*(u_n[i] - v_n[i]))


I guess that I am doing somewhere wrong while calculating the inner mesh points. Can anyone lemme know whats going wrong in here and ideas to rectify it?

• What is the stability criterion for selecting your timestep for this differential equation? Jan 21, 2019 at 20:58
• my stability criteria in this case is: F < = 0.5 Jan 21, 2019 at 21:23
• Double checking your equations, they seem to be non-linear but your FD looks like a linear iteration. Jan 21, 2019 at 21:41
• My advice, use FD to transform your PDEs to a system of transient ODEs and use scipy's odeint to integrate in time. It is optimized and robust, it will be difficult to optimize your own time-marching scheme to do better. Why reinvent the wheel? Jan 22, 2019 at 19:41
• @nluigi Do you have an example or some implementations for such coupled PDEs? Jan 23, 2019 at 17:27

This doesn't answer your question directly but instead suggests an alternate strategy. In general, unless you are interested in experimenting with numerical methods, I recommend using an existing PDE solver instead of trying to write one yourself. This advice is particularly true for numerically-challenging problems, like yours appears to be.

Since you are interested in Python, I have written a FiPy version of your problem, as shown below.

Unfortunately, I can't say that FiPy easily obtained a solution. But I have done some numerical experimentation and made a few observations. First, I should acknowledge that I don't understand much about the physics of these equations.

1. I changed your boundary conditions to zero-flux at the ends to be consistent with your initial conditions.

2. The laser pulse doesn't seem to be adding much heat to the system; the temperatures as a function of time change very little from the initial values. The temperature changes are so small, I have been plotting the difference in temperature compared with the initial conditions.

3. The large value of Gel appears to make the problem very numerically stiff and more difficult to solve. I've been experimenting with values several orders of magnitude less.

The Python code using the FiPy PDE solver is:

import numpy as np
import fipy
# sample parameter Gold
la0 = 429 # conductivity in W/mK
gma = 62.8 # thermal constant J/m^3K^2
Cl = 2.6e6 # phonon specific heat  in J/m^3K

Gel = 3.5e16 # lelectron phonon coupling constant
#Gel=1e10 # reduce value to improve numerics
tau_e = 0.04 # e relaxation t const
tau_l = 0.6 # l relaxation t const

Tm = 300.
L = 100 # sample thickness in nm
# Laser pulse parameters
T = 6.
T=1

R = 0.93 # reflectivity
I0 = 100 # fulence in J/m^2
tp = 0.1 # in ps
zz0 = 15.3 # penetration depth in nm

#F = 0.5

def las(x,t):
pt = 0.93*((1-R)/(zz0*tp))*I0*np.exp(-(x/zz0)-2.772*((t-2*tp)**2/tp**2))
return pt

def I(x):
tin = Tm
return tin

a = la0 # Ked

nx = 100
dx = L/nx

mesh = fipy.Grid1D(nx = nx, dx = dx)
x = mesh.cellCenters
#pt25 = mesh._getNearestCellID(([25.],))
u = fipy.CellVariable(name="u", mesh=mesh)
v = fipy.CellVariable(name="v", mesh=mesh)
valueLeft = 0
valueRight = 0
applyDirichletConstraint=False
if(applyDirichletConstraint):
u.constrain(valueRight, mesh.facesRight)
u.constrain(valueLeft, mesh.facesLeft)
v.constrain(valueRight, mesh.facesRight)
v.constrain(valueLeft, mesh.facesLeft)
else:
pi=np.math.pi
#Gumv = Gel*(u-v) # explicit source terms are not appropriate for this stiff system
Gumv = fipy.ImplicitSourceTerm(Gel, var=u) - fipy.ImplicitSourceTerm(Gel, var=v)
time = fipy.Variable()

eqn0 = fipy.TransientTerm(gma*u, var=u) == fipy.DiffusionTerm(a, var=u) - Gumv + las(x,time)
eqn1 = fipy.TransientTerm(Cl, var=v) == fipy.DiffusionTerm(a, var=v) + Gumv
eqn = eqn0 & eqn1
steps = 5000
dt = T/steps

ival = I(x)

u.setValue(ival)
v.setValue(ival)
centerX=([L/2.],)
uCenter = [u(centerX)]
vCenter = [v(centerX)]
ts  = [0.]
t = 0
for step in range(steps):
time.setValue(t+dt/2.)
eqn.solve(dt=dt)
uCenter.append(u(centerX))
vCenter.append(v(centerX))
t += dt
ts.append(t)

# plot difference from intial value
uCenter = np.asarray(uCenter) - Tm
vCenter = np.asarray(vCenter) - Tm

import matplotlib.pyplot as plot
plot.subplot(211);
plot.title('nx=%d, dt=%12.3e, Gel=%g' % (nx, dt, Gel))
p1=plot.plot(ts, uCenter)
plot.grid(b='on')
plot.xlabel('Time')
plot.ylabel('u at center')

plot.subplot(212);
p2=plot.plot(ts, vCenter)
plot.grid(b='on')
plot.xlabel('Time')
plot.ylabel('v at center')

plot.savefig('plot.pdf')

plot.show()