Problem after discretizing a PDE to ODE and solve with DifferentialEquations.jl

Hi, I was reading https://www.stochasticlifestyle.com/solving-systems-stochastic-pdes-using-gpus-julia/

And I would like to know if it is possible to discretize a 4th derivative PDE to ODE: ddt(u) = ddddx(u) and solve with DifferentialEquations.jl.
I am trying to discretize the RHS of the equation and solve it. But I got an error with Instability detected. Aborting. I increased the mesh size and changed to different solvers, but it is still not working.

Any help would be appreciated!

using DiffEqOperators
using BandedMatrices
using DifferentialEquations

m = 0.5
L = pi/m
n = 200 # 0,1,…,n-1,n
N = n-2 #2,3,…,n-3,n-2
# mesh from 0 to n
xl = LinRange(0.0,L,n)
dx = xl[2] - xl[1]

A = BandedMatrix(
-2 => ones(N-2)(1),
-1 => ones(N-1)
(-4),
0 => ones(N)(6),
1 => ones(N-1)
(-4),
2 => ones(N-2)*(1)
)

A[1,1] = 1
A[1,2] = -4/5
A[1,3] = 1/5
A[end,end] = 1
A[end,end-1] = -4/5
A[end,end-2] = 1/5

function myode(ddu, du, u, p, t)
ddu[:] = A*u./(dx^4)
end

x0 = LinRange(0, L, n)
u0 = sin.(m .*x0)[2:end-1]
du0 = cos.(m .*x0)[2:end-1]
tmax = 1.
p = 0
tspan = (0,tmax)
prob = SecondOrderODEProblem(myode,du0,u0,tspan,p)
sol = solve(prob, RadauIIA5(), saveat=0:0.01:tmax)

Are you sure that’s a stable discretization of the PDE? Do a Von Neumann analysis.

For reference, look through the following:

1 Like