Why does pdepe throw an error depending on boundary condition parameters?

6 views (last 30 days)
I want to simulate 1D wave propagation by solving the (acoustic) wave equation using pdepe.
Here the mathematical problem:
And here the corresponding code:
%% Define PDE system
x = linspace(0,10,50); %descretized space
t = linspace(0,0.1,50); %descretized time
m = 0;
sol = pdepe(m,@homWaveEq,@homWaveEq_ic,@homWaveEq_bc,x,t);
p = sol(:,:,1); % acoustic pressure
%% plot time evolution
xlabel('Distance x')
ylabel('Time t')
%% functions
function [c,f,s] = homWaveEq(x,t,u,dudx) % Equation to solve
K = 142000; % bulk modulus air
rho = 1.2; % density air
c = [1; 1];
f = [0;
s = [u(2);
% ---------------------------------------------
function u0 = homWaveEq_ic(x) % Initial Conditions
u0 = [0; 0];
% ---------------------------------------------
function [pl,ql,pr,qr] = homWaveEq_bc(xl,ul,xr,ur,t) % Boundary Conditions
freq = 40;
pl = [ul(1)-sin(2*pi()*freq*t); ul(2)-2*pi()*freq*cos(2*pi()*freq*t)];
ql = [0; 0];
pr = [0; 0];
qr = [1; 1];
The code works well for freq <= 40.
However, when freq > 40, I get the following error message:
"Error using pdepe (line 293)
Spatial discretization has failed. Discretization supports only parabolic and elliptic equations, with flux
term involving spatial derivative."
Has anyone an idea what the reason is for that and how to avoid it? I assume the error lies in the ode15s solver.
  1 Comment
Torsten on 19 Jan 2022
Edited: Torsten on 19 Jan 2022
I'm surprised that pdepe can solve this equation at all. Usually, the boundary condition to the right for u1 (which gives 0=0, thus something undetermined) leads to numerical problems.
Your equation is hyperbolic in nature and pdepe is designed to solve parabolic-elliptic problems.
You should use a solver for hyperbolic problems, e.g. CLAWPACK available at

Sign in to comment.

Accepted Answer

Bill Greene
Bill Greene on 20 Jan 2022
The ODE solver(ode15s) is numerically approximating the
iteration matrix for the PDE and, due to roundoff
error, deciding this is an ODE system that it cannot solve.
One way to get around this problem is to insure that your
initiation and boundary conditions are consistent. This is
often not strictly necessary, as you have noticed, but is
usually a good idea.
If you define freq in your main program and replace your
initial condition function with
function u0 = homWaveEq_ic(x) % Initial Conditions
if x==0
u0=[0; -2*pi*freq];
u0 = [0; 0];
you can obtain solutions for larger driving frequencies.

More Answers (0)




Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!