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