Appropriate method for solving coupled pdes

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.

Answer 1

[Bill Greene]

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

[Issy]

Hello Bill. Thank you very much for your answer and for the code. This will save me considerable time. It works excellently and i can easily adapt it to my needs. It is surprising that Mathworks describes pdepe as it does and i hope your comment about working for some hyperbolic pdes is useful to others. Many thanks once again and best wishes for the New Year. Ian.

[Issy]

Just a quick one please Bill. Should f=[0 0 a^2*dUdx(1) a^2*b^2*dUdx(2)]' ? Many thanks

[Bill Greene]

Oh yes, sorry, thanks for catching that! I corrected the code in my post above.

[Issy]

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]

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.