Why does pdepe throw an error depending on boundary condition parameters?
5 次查看(过去 30 天)
显示 更早的评论
Hello,
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
figure
surf(x,t,p)
xlabel('Distance x')
ylabel('Time t')
title('p(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;
K/rho*dudx(1)];
s = [u(2);
0];
end
% ---------------------------------------------
function u0 = homWaveEq_ic(x) % Initial Conditions
u0 = [0; 0];
end
% ---------------------------------------------
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];
end
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 个评论
Torsten
2022-1-19
编辑:Torsten
2022-1-19
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
采纳的回答
Bill Greene
2022-1-20
Interesting.
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];
else
u0 = [0; 0];
end
end
you can obtain solutions for larger driving frequencies.
更多回答(0 个)
另请参阅
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!