# Using OrdinaryDiffEq.jl solving 1d Wave Equation

**URL:** <https://discourse.julialang.org/t/using-ordinarydiffeq-jl-solving-1d-wave-equation/73782>\
**Category:** New to Julia\
**Tags:** ordinarydiffeq\
**Created:** [December 29, 2021, 7:44pm UTC](https://discourse.julialang.org/t/using-ordinarydiffeq-jl-solving-1d-wave-equation/73782 "2021-12-29T19:44:13Z")\
**Posts on this page:** 17\
**Page:** 1

<div class="post-metadata">

**Author:** ![CollinsPeriodLee](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/collinsperiodlee/32/23211_2.png) [@CollinsPeriodLee](https://discourse.julialang.org/u/CollinsPeriodLee)\
**Post date:** [December 29, 2021, 7:44pm UTC](https://discourse.julialang.org/t/using-ordinarydiffeq-jl-solving-1d-wave-equation/73782/1 "2021-12-29T19:44:13Z")

</div>

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

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

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

```

Any suggestion will be appreciated.

---

<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:** [December 29, 2021, 7:50pm UTC](https://discourse.julialang.org/t/using-ordinarydiffeq-jl-solving-1d-wave-equation/73782/2 "2021-12-29T19:50:51Z")

</div>

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](https://mitmath.github.io/18337/lecture14/pdes_and_convolutions). It’s a very classical result. I would try something like an upwinding discretization in space and see how that turns out.

---

<div class="post-metadata">

**Author:** ![CollinsPeriodLee](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/collinsperiodlee/32/23211_2.png) [@CollinsPeriodLee](https://discourse.julialang.org/u/CollinsPeriodLee)\
**Post date:** [December 29, 2021, 7:55pm UTC](https://discourse.julialang.org/t/using-ordinarydiffeq-jl-solving-1d-wave-equation/73782/3 "2021-12-29T19:55:58Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![CollinsPeriodLee](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/collinsperiodlee/32/23211_2.png) [@CollinsPeriodLee](https://discourse.julialang.org/u/CollinsPeriodLee)\
**Post date:** [December 29, 2021, 8:31pm UTC](https://discourse.julialang.org/t/using-ordinarydiffeq-jl-solving-1d-wave-equation/73782/4 "2021-12-29T20:31:26Z")

</div>

If I remove the dot,

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

```

It works… Is there any possible explanation?

---

<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:** [December 29, 2021, 8:40pm UTC](https://discourse.julialang.org/t/using-ordinarydiffeq-jl-solving-1d-wave-equation/73782/5 "2021-12-29T20:40:51Z")

</div>

> [@CollinsPeriodLee](#):
>
> 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`.

---

<div class="post-metadata">

**Author:** ![CollinsPeriodLee](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/collinsperiodlee/32/23211_2.png) [@CollinsPeriodLee](https://discourse.julialang.org/u/CollinsPeriodLee)\
**Post date:** [December 29, 2021, 10:51pm UTC](https://discourse.julialang.org/t/using-ordinarydiffeq-jl-solving-1d-wave-equation/73782/6 "2021-12-29T22:51:01Z")

</div>

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

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

```

---

<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:** [December 29, 2021, 10:51pm UTC](https://discourse.julialang.org/t/using-ordinarydiffeq-jl-solving-1d-wave-equation/73782/7 "2021-12-29T22:51:49Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![CollinsPeriodLee](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/collinsperiodlee/32/23211_2.png) [@CollinsPeriodLee](https://discourse.julialang.org/u/CollinsPeriodLee)\
**Post date:** [December 29, 2021, 11:58pm UTC](https://discourse.julialang.org/t/using-ordinarydiffeq-jl-solving-1d-wave-equation/73782/8 "2021-12-29T23:58:22Z")

</div>

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

---

<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:** [December 30, 2021, 2:17am UTC](https://discourse.julialang.org/t/using-ordinarydiffeq-jl-solving-1d-wave-equation/73782/9 "2021-12-30T02:17:24Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![ranocha](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ranocha/32/35588_2.png) [@ranocha](https://discourse.julialang.org/u/ranocha)\
**Post date:** [December 30, 2021, 7:33am UTC](https://discourse.julialang.org/t/using-ordinarydiffeq-jl-solving-1d-wave-equation/73782/10 "2021-12-30T07:33:37Z")

</div>

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](https://ranocha.de/SummationByPartsOperators.jl/stable/tutorials/wave_equation).

---

<div class="post-metadata">

**Author:** ![CollinsPeriodLee](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/collinsperiodlee/32/23211_2.png) [@CollinsPeriodLee](https://discourse.julialang.org/u/CollinsPeriodLee)\
**Post date:** [December 30, 2021, 4:39pm UTC](https://discourse.julialang.org/t/using-ordinarydiffeq-jl-solving-1d-wave-equation/73782/11 "2021-12-30T16:39:53Z")

</div>

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

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

---

<div class="post-metadata">

**Author:** ![CollinsPeriodLee](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/collinsperiodlee/32/23211_2.png) [@CollinsPeriodLee](https://discourse.julialang.org/u/CollinsPeriodLee)\
**Post date:** [December 30, 2021, 4:40pm UTC](https://discourse.julialang.org/t/using-ordinarydiffeq-jl-solving-1d-wave-equation/73782/12 "2021-12-30T16:40:34Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![CollinsPeriodLee](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/collinsperiodlee/32/23211_2.png) [@CollinsPeriodLee](https://discourse.julialang.org/u/CollinsPeriodLee)\
**Post date:** [December 30, 2021, 4:46pm UTC](https://discourse.julialang.org/t/using-ordinarydiffeq-jl-solving-1d-wave-equation/73782/13 "2021-12-30T16:46:14Z")

</div>

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 ?

---

<div class="post-metadata">

**Author:** ![ranocha](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ranocha/32/35588_2.png) [@ranocha](https://discourse.julialang.org/u/ranocha)\
**Post date:** [December 30, 2021, 6:26pm UTC](https://discourse.julialang.org/t/using-ordinarydiffeq-jl-solving-1d-wave-equation/73782/14 "2021-12-30T18:26:54Z")

</div>

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](https://ranocha.de/SummationByPartsOperators.jl/stable/api_reference/#SummationByPartsOperators.Mattsson2012)  
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.

---

<div class="post-metadata">

**Author:** ![CollinsPeriodLee](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/collinsperiodlee/32/23211_2.png) [@CollinsPeriodLee](https://discourse.julialang.org/u/CollinsPeriodLee)\
**Post date:** [December 30, 2021, 7:59pm UTC](https://discourse.julialang.org/t/using-ordinarydiffeq-jl-solving-1d-wave-equation/73782/15 "2021-12-30T19:59:34Z")

</div>

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)](https://ranocha.de/SummationByPartsOperators.jl/stable/tutorials/variable_linear_advection/), right? I also need to implement the stable boundary procedures.

---

<div class="post-metadata">

**Author:** ![ranocha](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ranocha/32/35588_2.png) [@ranocha](https://discourse.julialang.org/u/ranocha)\
**Post date:** [December 31, 2021, 7:12am UTC](https://discourse.julialang.org/t/using-ordinarydiffeq-jl-solving-1d-wave-equation/73782/16 "2021-12-31T07:12:35Z")

</div>

That’s basically it. You can create a variable-coefficient derivative operator using [`var_coef_derivative_operator`](https://ranocha.de/SummationByPartsOperators.jl/stable/api_reference/#SummationByPartsOperators.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.

---

<div class="post-metadata">

**Author:** ![CollinsPeriodLee](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/collinsperiodlee/32/23211_2.png) [@CollinsPeriodLee](https://discourse.julialang.org/u/CollinsPeriodLee)\
**Post date:** [December 31, 2021, 4:32pm UTC](https://discourse.julialang.org/t/using-ordinarydiffeq-jl-solving-1d-wave-equation/73782/17 "2021-12-31T16:32:15Z")

</div>

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