# Using OrdinaryDiffEq.jl solving 1d Wave Equation

Hi, I’m using OrdinaryDiffEq.jl to solve the wave Equation, but there is Instability detected. Attached is my code

``````using OrdinaryDiffEq, Sundials
using LinearAlgebra
using Plots; gr()

function vright(N,ϵ, xs,dx)
d = -2 * ones(N-2)  # main diagonal
du = ones(N - 3) # off diagonal
D = (ϵ^2/dx^2) * diagm(-1 => du, 0 => d, 1 => du)
function (du,u,p,t)
mul!(du,D,u)
du .= du .+ rhs.(xs)
end
end

function rhs(x)
return exp(x) - exp(-x)
end

# Construct the problem
function waves(N, ϵ, xrange)
dx = (xrange[2] - xrange[1])/(N-1)
xs = (1:N-2)*dx .+ xrange[1]

f1= vright(N,ϵ, xs, dx)

μ0 = 0.0; σ0 = 0.25
f0 = x -> x #((x-μ0)/σ0)^2
f0prime = x -> sin(x)
u0 = f0.(xs)
du0 = f0prime.(xs)

prob = SecondOrderODEProblem(f1,du0, u0, (0.0, 1.0))
return xs, prob
end;

function exact(x,t)
return x + 1/3*sin(x)*sin(3t) -2/9*sinh(x)+2/9*sinh(x)*cosh(3t)
end;

xrange = [-1,1]
xs, prob = waves(1001, 3.0,xrange);
sol = solve(prob, DPRKN12(); abstol=1e-14, reltol=1e-14);
``````

The error reported is

``````┌ Warning: Instability detected. Aborting
└ @ SciMLBase ~/.julia/packages/SciMLBase/mndcy/src/integrator_interface.jl:351
``````

Any suggestion will be appreciated.

1 Like

Are you should that choice of spatial discretization is stable. For the advection equation I know you can prove that the choice you made here is unconditionally unstable using a Von Neumann analysis similar to the one shown in PDEs, Convolutions, and the Mathematics of Locality. It’s a very classical result. I would try something like an upwinding discretization in space and see how that turns out.

Got it. Thanks, I will have a try and update!

If I remove the dot,

``````    function (du,u,p,t)
du = D*u .+ rhs.(xs)
return du
end
``````

It works… Is there any possible explanation?

Yes, because then the derivative is just zero because it’s not mutating. But then the solution is constant, which isn’t what you want. You should probably just make that `return nothing`.

1 Like

Hi, Chris, I used the upwind scheme, but the warning that Instability detected still popped out.

``````function vright(N,ϵ, xs,dx)
d = -2 * ones(N-2)  # main diagonal
du = ones(N - 3) # off diagonal
A = UpwindDifference{Float64}(2,2,dx,N-2,-1)
function (du, u,p,t)
bc = DirichletBC(0.,0.)
mul!(du,A,bc*u)
@. du = du + rhs(xs)
end
end

function rhs(x)
return exp.(x) .- exp.(-x)
end

# Construct the problem
function waves(N, ϵ, xrange)
dx = (xrange[2] - xrange[1])/(N-1)
xs = (1:N-2)*dx .+ xrange[1]

f1= vright(N,ϵ, xs, dx)

f0 = x -> x #((x-μ0)/σ0)^2
f0prime = x -> sin(x)
u0 = f0.(xs)
du0 = f0prime.(xs)

prob = SecondOrderODEProblem(f1,du0, u0, (0.0, 1.0))
return xs, prob
end;

function exact(x,t)
return x + 1/3*sin(x)*sin(3t) -2/9*sinh(x)+2/9*sinh(x)*cosh(3t)
end;

xrange = [-10,10]
xs, prob = waves(2001, 3.0,xrange);
sol = solve(prob, DPRKN12(); abstol=1e-14, reltol=1e-14);
``````

What about with other ODE solvers like Vern9? Did you try any methods for stiff equations like QNDF?

I replace the ODE solvers to Vern9, but the result is incorrect, I’m checking that.

Just write out the upwinding by hand. DiffEqOperators isn’t ready for prime time.

What kind of scheme and which boundary conditions do you want to use? You could look at my package SummationByPartsOperators.jl - I have implemented some stable finite difference methods there, see the brief tutorial on the 1D wave equation.

4 Likes

Hi, Chris, Thanks for the advice, I just figured the issue out. It is the input of the function f1.

``````function vright(N,ϵ, xs,dx)
d = -2 * ones(N-2)  # main diagonal
du = ones(N - 3) # off diagonal
D = (ϵ^2/dx^2) * diagm(-1 => du, 0 => d, 1 => du)
D
end

function f1!(ddu, du, u,p,t)
D, xs = p
bc = DirichletBC(0.,0.)
ddu .= D * bc*u .+ rhs(xs)
end
``````

Generally, for `SecondOrderODEProblem{isinplace}(f,du0,u0,tspan,callback=CallbackSet())`, the function `f` needs either 4 inputs `du,u,p,t` or 5 inputs `ddu,du,u,p,t` for the in-place case. In the previous case, I use in-place 4 inputs that causes `Instability Detected`.

On the other hand, I got O(1e-7) absolute error when `t=0.1` and O(1e-3) absolute error at `t=1`. Is this what the unconditionally instability you mentioned?

Got it! Thanks, I will have a look and give the feedback here!

I just tested the package, it is just what I really want! The boundary conditions is Dirichlet boundary condition. I mainly focus on getting the solutions of wave equation, thus I could use any schemes.

I have a further question, is it possible solving the wave equation with variable coefficients ?

Yes, that’s possible. You basically need a second-derivative operator with variable coefficients - I have implemented such a source of coefficients, see API reference · SummationByPartsOperators.jl
If I remember correctly, you can’t use the code from the tutorial linked above directly. However, Ken Mattsson et al. should have a paper describing stable boundary procedures for the case you are interested in. You should find some references in the source code implementing the methods shown in the tutorial.

1 Like

Thanks for your kind explanation. Please correct me if I miss any part.

I need to implement a semidiscretization like `SummationByPartsOperators.VariableLinearAdvectionNonperiodicSemidiscretization`, and call it as in Linear advection equation with variable coefficients · SummationByPartsOperators.jl (ranocha.de), right? I also need to implement the stable boundary procedures.

That’s basically it. You can create a variable-coefficient derivative operator using `var_coef_derivative_operator`. Note that you need a second-derivative operator as in the wave equation tutorial linked above, not a first-derivative operator as in the linear advection tutorial.
You are of course welcome to make a PR to add the wave equation with variable coefficients as tutorial to the documentation - and get some feedback there.

1 Like

OK, once I finish, I will make a PR. Thank you very much!