# Numerical instability of partial differential equations

**URL:** <https://discourse.julialang.org/t/numerical-instability-of-partial-differential-equations/39603>\
**Category:** General Usage\
**Tags:** diffeq\
**Created:** [May 16, 2020, 10:24pm UTC](https://discourse.julialang.org/t/numerical-instability-of-partial-differential-equations/39603 "2020-05-16T22:24:39Z")\
**Posts on this page:** 7\
**Page:** 1

<div class="post-metadata">

**Author:** ![Chong\_Wang](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chong_wang/32/20307_2.png) [@Chong\_Wang](https://discourse.julialang.org/u/Chong_Wang)\
**Post date:** [May 16, 2020, 10:24pm UTC](https://discourse.julialang.org/t/numerical-instability-of-partial-differential-equations/39603/1 "2020-05-16T22:24:39Z")

</div>

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

```julia
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

```julia
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[1] = (ys[2]-ys[1])/(dx)
    dys[end] = (ys[end]-ys[end-1])/(dx)
    return dys
end

```

The evolution of the state is implemented as

```julia
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

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

```

I obtain something like this

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

```

 ![image](https://global.discourse-cdn.com/julialang/original/3X/3/2/322c10c7897c36d57ab00482b6cf9f4b7b4a9e93.png)  
The solution is oscillating like crazy in some region.

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

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

```

 ![image](https://global.discourse-cdn.com/julialang/original/3X/e/3/e32e02d046233d068637e9b3530d37ab5fbbc63e.png)

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](https://en.wikipedia.org/wiki/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.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [May 17, 2020, 1:26am UTC](https://discourse.julialang.org/t/numerical-instability-of-partial-differential-equations/39603/2 "2020-05-17T01:26:46Z")

</div>

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.

> [@Chong\_Wang](#):
>
> dys[i] = (ys[i+1]-ys[i-1])/(2dx)

This is a second order approximation to the forward derivative.

> [@Chong\_Wang](#):
>
> differentiate(differentiate(ms, dx), dx)

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 `dt`s 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)

---

<div class="post-metadata">

**Author:** ![Chong\_Wang](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chong_wang/32/20307_2.png) [@Chong\_Wang](https://discourse.julialang.org/u/Chong_Wang)\
**Post date:** [May 17, 2020, 3:17am UTC](https://discourse.julialang.org/t/numerical-instability-of-partial-differential-equations/39603/3 "2020-05-17T03:17:32Z")

</div>

Thank you for the reply.

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)

```julia
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()

```

 ![image](https://global.discourse-cdn.com/julialang/original/3X/2/b/2bf1ae5597b5e5bfd608957ee6bee0ff435b0947.png)

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

```julia
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

```julia
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]

```

 ![true](https://global.discourse-cdn.com/julialang/original/3X/5/e/5e74506c62c73f0e8ed7fe2d656e2ec9e10bab76.png)  
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.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [May 17, 2020, 4:02am UTC](https://discourse.julialang.org/t/numerical-instability-of-partial-differential-equations/39603/4 "2020-05-17T04:02:30Z")

</div>

> [@Chong\_Wang](#):
>
> Also, reducing `dt` doesn’t work; actually, the result is worse for

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

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

---

<div class="post-metadata">

**Author:** ![YingboMa](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/yingboma/32/2181_2.png) [@YingboMa](https://discourse.julialang.org/u/YingboMa)\
**Post date:** [May 17, 2020, 4:11am UTC](https://discourse.julialang.org/t/numerical-instability-of-partial-differential-equations/39603/5 "2020-05-17T04:11:10Z")

</div>

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.

```julia
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:

 ![plt](https://global.discourse-cdn.com/julialang/original/3X/9/8/98aa3a57d63c69d887977fb4407d141593af532d.png)

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

---

<div class="post-metadata">

**Author:** ![Chong\_Wang](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chong_wang/32/20307_2.png) [@Chong\_Wang](https://discourse.julialang.org/u/Chong_Wang)\
**Post date:** [May 17, 2020, 5:03am UTC](https://discourse.julialang.org/t/numerical-instability-of-partial-differential-equations/39603/6 "2020-05-17T05:03:50Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![YingboMa](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/yingboma/32/2181_2.png) [@YingboMa](https://discourse.julialang.org/u/YingboMa)\
**Post date:** [May 17, 2020, 5:15am UTC](https://discourse.julialang.org/t/numerical-instability-of-partial-differential-equations/39603/7 "2020-05-17T05:15:40Z")

</div>

> [@Chong\_Wang](#):
>
> 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.

`Tsit5` has a larger stability region than the Euler method.

> [@Chong\_Wang](#):
>
> Actually, I tried `DifferentialEquations.jl` , but to conserve the length of m\boldsymbol{m} at every spatial point, I used a callback function and it takes forever to finish.

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