# Numerical instability of partial differential equations

Probably I should have asked this on a mathematical forum, but since I will present julia code and there are a lot of experts in this forum, probably I can get better help here.

I am trying to solve a differential equation

\partial_t \boldsymbol{m}=\partial_x^2\boldsymbol{m}-(\boldsymbol{m}\cdot\partial_x^2\boldsymbol{m})\boldsymbol{m},

where \boldsymbol{m} is a unit vector, which is clearly preserved by the equation.

I start from initial state \phi_0(x)=0 and \theta_0(x)=2\arctan \exp (x/5), where \theta_0 and \phi_0 are spherical coordinates:

\boldsymbol{m}_0(x)=(\sin \theta_0 \cos \phi_0,\sin \theta_0 \sin \phi_0,\cos \theta_0).

The corresponding code is

function construct_m(θs, ϕs)
ms = zeros(3, length(θs))
for i in 1:length(θs)
θ = θs[i]; ϕ = ϕs[i]
ms[:, i] += [sin(θ)cos(ϕ), sin(θ)sin(ϕ), cos(θ)]
end
return ms
end

xs = -250:0.01:250
ms = construct_m(2*atan.(exp.(xs./5)), zeros(length(xs)));


where I am sampling on [-250, 250] with 0.01 accuracy.

The numerical differentiation is implemented as

function differentiate(ys, dx)
dys = zeros(size(ys))
for i in 2:length(ys)-1
dys[i] = (ys[i+1]-ys[i-1])/(2dx)
end
dys = (ys-ys)/(dx)
dys[end] = (ys[end]-ys[end-1])/(dx)
return dys
end


The evolution of the state is implemented as

using LinearAlgebra
function get_next_m(ms, dx, dt=0.1)
N = size(ms, 2)
dms = zeros(3, N)
for α in 1:3
dms[α, :] = differentiate(differentiate(ms[α, :], dx), dx)
end
nms = zeros(3, N)
for i in 1:N
m = ms[:, i]; dm = dms[:, i]
nms[:, i] = normalize(m+(dm-m*(m⋅dm))*dt)
end
return nms
end


However, only after 8 steps

for i in 1:8
ms = get_next_m(ms, 0.01, 0.01)
end


I obtain something like this

using PyPlot
figure(figsize=(15, 3))
plot(xs, ms[3, :])
xlim(-200, 200)


The solution is oscillating like crazy in some region.

However, if I only run for 7 steps, it still looks fine

for i in 1:7
ms = get_next_m(ms, 0.01, 0.01)
end


Thoughts:

• I believe this instability is due to numerical error in finite difference method. By adding the derivative \partial_x \boldsymbol{m} to \boldsymbol{m} and then differentiate \boldsymbol{m} again in the next iteration, I am differentiating \boldsymbol{m} up to 2500 order, which is clearly numerical unstable for finite difference.
• However, I believe I am following the standard procedure in solving such an equation. I scanned through a few online tutorials and everyone seems to use finite difference method for similar equations
• I noticed some type of differential equations are stiff equation. However, according to the wikipedia, stiff equations are oscillating with respect to t. Therefore I am not sure this is related.

I am wondering whether my equation is special in a way it is more numerical unstable than, for example, heat equation. And how can I solve this instability for my solution?

Edit:
This is actually a MWE for my actual code. Unfortunately, I made a mistake in reducing my code to MWE, and I have updated the code and results in this revision. After fixing the mistake, this MWE behaves more like my actual code: just after several steps, it blows off unexpectedly.

It’s a bit hard to analyze the stability quickly here, so I’ll at least jot down a few things that comes to mind to look at and try.

This is a second order approximation to the forward derivative.

You then do it twice. What you get is a valid scheme, but it’s not necessarily the one most people would try on this kind of thing. It’s not symmetric so it’s probably less stable then the standard scheme when applied to even just the heat equation part. Try using the standard u[i-1] - 2u[i] + u[1+1] central difference on both terms and see how that terms out. The second term might do better with a kind of upwinding-ish thing that you have going on here, so you might want to try 2nd order central on the first term and your current on the second term.

If that doesn’t work out, try using smaller dts and, if that doesn’t work out, see what’s happening in 10 steps. You should be able to figure out why terms in the middle there are continuing to grow unstably just by following the steps of the solution. Since it’s not on all of the downward portion and only specifically past zero, that’s quite an interesting hint. Off the top of my head, my guess would be that the instability of the equation stems from the fact that your update equation at u[i] doesn’t have u[i] in it, so you get an instability similar to leapfrog methods where you have terms oscillating between even and odd points but never hitting the middle, and this oscillation is then just growing (which then points to using a different discretization like the one I mentioned)

2 Likes

Unfortunately, I made a mistake in my original code. I have updated the post to reflect this.

After using the central difference approximation for second order derivative, the behavior is actually a little bit worse.

The full code (with central difference)

using LinearAlgebra, PyPlot
function differentiate_twice(ys, dx)
dys = zeros(size(ys))
for i in 2:length(ys)-1
dys[i] = (ys[i+1]+ys[i-1]-2ys[i])/(dx^2)
end
return dys
end
function construct_m(θs, ϕs)
ms = zeros(3, length(θs))
for i in 1:length(θs)
θ = θs[i]; ϕ = ϕs[i]
ms[:, i] += [sin(θ)cos(ϕ), sin(θ)sin(ϕ), cos(θ)]
end
return ms
end
function get_next_m(ms, dx, dt)
N = size(ms, 2)
dms = zeros(3, N)
for α in 1:3
dms[α, :] = differentiate_twice(ms[α, :], dx)
end
nms = zeros(3, N)
for i in 1:N
m = ms[:, i]; dm = dms[:, i]
nms[:, i] = normalize(m+(dm-m*(m⋅dm))*dt)
end
return nms
end
xs = -250:0.01:250
ms = construct_m(2*atan.(exp.(xs./5)), zeros(length(xs)))
for i in 1:6
ms = get_next_m(ms, 0.01, 0.01)
end
figure(figsize=(15, 3))
plot(xs, ms[3, :])
xlim(-100, 100)
tight_layout()


I tried to trace the program but cannot figure out why, partly because the oscillation appears so suddenly that it is difficult to find a hint. At the 5th step, \boldsymbol{m} still looks fine.

Also, reducing dt doesn’t work; actually, the result is worse for

for i in 1:60
ms = get_next_m(ms, 0.01, 0.001)
end


BTW, Mathematica correctly solves this problem by the following (ugly) code

sol = NDSolve[{D[mx[t, x], t] ==
D[mx[t, x], {x, 2}] -
mx[t, x] (D[mx[t, x], {x, 2}] mx[t, x] +
D[mz[t, x], {x, 2}] mz[t, x]),
D[mz[t, x], t] ==
D[mz[t, x], {x, 2}] -
mz[t, x] (D[mx[t, x], {x, 2}] mx[t, x] +
D[mz[t, x], {x, 2}] mz[t, x]),
mx[0, x] == Sin[2 ArcTan[Exp[x/5]]],
mz[0, x] == Cos[2 ArcTan[Exp[x/5]]],
mx[t, -250] == Sin[2 ArcTan[Exp[-250/5]]],
mx[t, 250] == Sin[2 ArcTan[Exp[250/5]]],
mz[t, -250] == Cos[2 ArcTan[Exp[-250/5]]],
mz[t, 250] == Cos[2 ArcTan[Exp[250/5]]]}, {mx, mz}, {t, 0,
250}, {x, -250, 250}, AccuracyGoal -> 6, PrecisionGoal -> 6,
StartingStepSize -> 0.05]


where the blue (yellow) line is m_z at t=0 (t=250). This equation describes an expansion of the domain wall. Unfortunately, not much insight can be gained since I don’t even know what method Mathematica is choosing for the problem.

That’s a sign of instability in the spatial discretization.

Maybe it makes more sense to do this in Fourier space?

1 Like

I modified your code a little to be more idiomatic Julia, and used the 5-th order Tsit5 integrator from OrdinaryDiffEq.jl to do the time-stepping. Also, I added some comments in the code.

using OrdinaryDiffEq, Plots
function differentiate_twice!(dys, ys, dx) # inplace operation saves memory
@inbounds for i in 2:length(ys)-1
dys[i] = (ys[i+1]+ys[i-1]-2ys[i])/(dx^2)
end
return dys
end
function rhs!(dms, ms, dx, t)
(size(dms, 1) == size(ms, 1) && size(dms, 2) == size(ms, 2) == 3) || throw(DimensionMismatch("Incorrect dimension of dms or ms."))
N = size(ms, 1)
# julia is column major, so iterating on the leading dimension is faster
for j in 1:3 # scalar index usually uses lowercase Latin letters
differentiate_twice!(@view(dms[:, j]), @view(ms[:, j]), dx)
end
# column major iteration
@inbounds for j in 1:3, i in 1:N
dm = dms[i, j]
m = dm[i, j]
# tell the compiler that is ok to fuse and reorder floating point calculations, but @fastmath should be used with care.
dms[i, j] = @fastmath dm - (m * dm) * m
end
return nothing
end
function construct_m(θs, ϕs)
ms = zeros(length(θs), 3)
@inbounds for i in 1:length(θs)
θ = θs[i]; ϕ = ϕs[i]
sθ, cθ = sincos(θ) # fused sin/cos calculation
sϕ, cϕ = sincos(ϕ)
ms[i, 1] += sθ*cϕ
ms[i, 2] += sθ*sϕ
ms[i, 3] += cθ
end
return ms
end
dx = 0.1 # 0.01 is too small, and it leads to stiffness
xs = -250:dx:250
ms = construct_m(2*atan.(exp.(xs./5)), zeros(length(xs)))
prob = ODEProblem(rhs!, ms, (0.0, 100.0), dx)
sol = solve(prob, Tsit5())
plt1 = plot(xs, sol(0.0), title="t=0")
plt2 = plot(xs, sol(25.0), title="t=25")
plt3 = plot(xs, sol(50.0), title="t=50")
plt4 = plot(xs, sol(100.0), title="t=100")
plot(plt1, plt2, plt3, plt4, dpi=200)
savefig("plt.png")


Plot:

Also, as Chris suggested, solving this PDE in Fourier space would be better.

4 Likes

Thank you.

I still don’t understand what is the problem. If it is the spatial discretization, I changed my dx from 0.01 to 0.1 (like what is done by @YingboMa ) but the solution is still unstable. It seems Tsit5() in DifferentialEquations.jl is necessary, but Tsit5(), as far as I can understand, deals with the time derivative, not the spatial one. Also, when would I generally expect problems with spatial descretization?

Actually, I tried DifferentialEquations.jl, but to conserve the length of \boldsymbol{m} at every spatial point, I used a callback function and it takes forever to finish. @YingboMa In your plots, it seems that the length of \boldsymbol{m} is not conserved, at least at x=0. I ran your code with atol=1.0e-10 and rtol=1.0e-10 and the result seems to be the same. (BTW, I wasn’t aware DifferentialEquations.jl works for matrix input, so my code calling DifferentialEquations.jl was ugly)

I agree that discretization in Fourier space might help. I’ll try that.

Tsit5 has a larger stability region than the Euler method.

Conservation, in general, can better be formulated as a differential-algebraic equation.

3 Likes