Appropriate method for solving coupled pdes

68 views (last 30 days)
Hello. I am looking to solve the following system of modified wave equations:
I have looked at pdepe and thought there might be an equivalent for hyperbolic equations. I thought of using the pde modeller in 2D with a thin domain and zero neumann boundary conditions on the sides but the equations did not fit the pde modeller (the single derivative time terms could not be accommodated in the standard form in pde modeller). Am I overlooking a method in matlab or simulink or do I have to write my own finite difference scheme? Many thanks for any pointers.

Accepted Answer

Bill Greene
Bill Greene on 31 Dec 2020
Edited: Bill Greene on 31 Dec 2020
Yes, although it is true that the documentation for pdepe describes it as a solver for parabolic systems, it can often obtain solutions to systems of hyperbolic PDE with no problems.
To convert your two-equation system to a form acceptable to pdepe, I defined two auxiliary variables to yield the following four-equation system
I made a few assumptions about the values of the constants, the initial, and boundary conditions and used the following MATLAB code to solve this system
function matlabAnswers_12_31_2020
nx=31;
x = linspace(0,1,nx);
nx2=ceil(nx/2);
t = linspace(0, 1, 50);
m = 0;
sol = pdepe(m,@pdefun,@pdeic,@pdebc,x,t);
u=sol(:,:,1);
v=sol(:,:,2);
figure; plot(x, u(end,:) ,x,v(end,:)); grid; legend('u', 'v');
title 'solution at final time';
figure; plot(t, u(:,nx2),t,v(:,nx2)); grid; legend('u', 'v');
title 'solution at the center as a function of time';
end
function [ c, f, s ] = pdefun ( x, t, U, dUdx)
a=1; b=1;
c=[1 1 1 1]';
f=[0 0 a^2*dUdx(1) a^2*b^2*dUdx(2)]';
s=[U(3) U(4) U(4) -U(3)]';
end
function u0 = pdeic(x)
u0=[sin(pi*x) sin(3*pi*x) 0 0]';
end
function [ pl, ql, pr, qr ] = pdebc ( xl, ul, xr, ur, t )
pl = ul;
ql = [0 0 0 0]';
pr = ur;
qr = [0 0 0 0]';
end
  5 Comments
Issy
Issy on 31 Dec 2020
Hi Bill. I have just checked the code against another method I used earlier ie. taking the fourier transform and solving analytically the first order system in the transformed space (this can be done in maple). I then used the 'integral' command in matlab to numerically invert the transform. The two results are the same which is good to know with all the numerics involved. Many thanks for your help once again. I hope your observation about solving hyperbolic pdes in matlab becomes more widely appreciated. Best wishes, Ian
Bill Greene
Bill Greene on 1 Jan 2021
Thanks for the check on the pdepe solution. When the boundary or initial conditions are discontinuous, solving hyperbolic PDE with pdepe can be challenging. I suspect this is why the documentation doesn't emphasize this class of problem.

Sign in to comment.

More Answers (0)

Tags

Products


Release

R2020a

Community Treasure Hunt

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

Start Hunting!