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.
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.
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 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.
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?
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.
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.