Solving coupled set of PDEs and ODEs for 1D-problem using pde1dM

Hello,

I'like to solve a coupled set of PDEs (N=4 unknowns) and ODEs (M=1 unknown) using the solver pde1dM. A detailed description of the mathematical problem is attached. The problem is an extension of the standard wave equation, which I already posted here and I was able to solve successfully with the pdepe solver.

The ODE should be solved at the same locations as the PDE. Therefore, I set the vector xOde (of length L) equal to meshPts. When I am running my code, the solver stops after T time steps because of not meeting the integration tolerances. The solution matrix of the PDE set (sol) has the size of TxLxN, but the solution matrix of the ODE (odeSol) is only of size Tx1xM but should be TxLxM.

Does anyone know the reason for that ?

Thanks in advance!

I am running the following code:

clear all

set(0,'DefaultLegendAutoUpdate','off')

set(0, 'DefaultTextInterpreter', 'latex')

set(0, 'DefaultLegendInterpreter', 'latex')

set(0, 'DefaultAxesTickLabelInterpreter', 'latex')

global freq

global exc

global par_A

global par_alpha

global par_B

global par_beta

global par_C

global par_lambda

global par_Kinf

global par_rhoInf

% addpath("pde1dM-master") %load pde-ode solver

% load('../RationalFunction_Fitting/MPP_Par_K-N1M0_rho-N0M1.mat') %load equivalent fluid model

par_A = 5.547518377073091e-04;

par_alpha = -2.545221986206693e+03;

par_B = -96.874082239056780;

par_beta = -2.946796763623190e+02;

par_C = -18.908895618063692;

par_lambda = -4.128054318161039e+03;

par_Kinf = 4.483030705842021e-07;

par_rhoInf = 0.037897808274425;

%% Input Parameters

freq = 100;

exc = 'sinburst' % excitation signal: 'harmonic' or 'sinburst'

%% Define PDE system

Nx = 50

meshPts = linspace(0,10,Nx);

t = linspace(0,0.1,200);

delta_t = t(2)-t(1);

xOde = meshPts; % common node positions of PDE and ODE (all of them are coupled)

%% solve pde

m = 0;

[sol,odeSol] = pde1dm(m,@pde_F,@pde_IC,@pde_BC,meshPts,t,@ode_F, @ode_IC,xOde);

p = sol(:,:,1);

%% plot solution

%time evolution

figure

for k = 1:length(t)

plot(meshPts,p(k,:))

% hold on

% plot (x,q(k,:)/(2*pi*40))

% hold off

xlabel('Distance x')

ylabel('Time t')

ylim([-3 3])

title('p(t) along axis')

%pause(0.2)

end

% 3d plot of solution

figure

surf(meshPts,t,p)

xlabel('Distance x')

ylabel('Time t')

title('p(t)')

%% functions

function [c,f,s] = pde_F(x,t,u,dudx,v,vDot) % PDE to be solved

global par_A

global par_alpha

global par_B

global par_beta

global par_C

global par_Kinf

global par_rhoInf

c = [1; 1; 1; 1];

f = [0;

par_rhoInf/par_Kinf*dudx(1)+ 1/par_Kinf*par_B*dudx(3) + 1/par_Kinf*par_C*dudx(4);

0;

0];

s = [u(2);

-par_A/par_Kinf*v(1);

-par_alpha*u(3)-par_beta*u(4)+u(1);

-par_alpha*u(4)+par_beta*u(3)];

end

% ---------------------------------------------

function u0 = pde_IC(x) % Initial conditions of PDE

if x==0

u0=[0; 0; 0; 0];

else

u0 = [0; 0; 0; 0];

end

end

% ---------------------------------------------

function [pl,ql,pr,qr] = pde_BC(xl,ul,xr,ur,t,v,vDot) % BCs

global freq

global exc

pl = [ul(1)-(1+cos(2*pi()*freq*t+pi));...

ul(2)+2*pi()*freq*sin(2*pi()*freq*t+pi);...

0;

0];

if t >= 1/freq & exc == 'sinburst' % otherwise harmonic excitation

pl = [ul(1); ul(2); 0; 0];

end

ql = [0; 0; 1; 1];

pr = [0; 0; 0; 0];

qr = [1; 1; 1; 1];

end

function R=ode_F(t,v,vdot,x,u,DuDx,f, dudt, du2dxdt) % ODE to be solved

global par_lambda

R= vdot +par_lambda*v(1)-dudt(2);

end

function value = ode_IC() % initial condition of ODE

value = 0;

end

Answer 1

[Bill Greene]

Because your equation (2) is a PDE and not an ODE (as I pointed out

in your previous post), it will be quite challenging to solve using

pde1dm treating it as a collection of ODE.

I strongly suggest you just add equation (2) as another PDE in your

system, i.e. a total of five PDE. For this PDE you can define

the flux equal to zero as boundary conditions. Because the flux is

equal to zero everywhere, this BC will not affect the solution but will

satisfy the requirements of pdepe (or pde1dm).

[Paul Maurerlehner]

For me equation (2) is no PDE, because no spatial derivatives occur, but

The problem of your suggestion is, that I would need the second time derivative of of the pressure for equation (2), which I don't have available in the PDE function. Order reduction by substituation as i am already doint is not possible, because the second time derivative (time derivative of u_2) of the pressure is already used for equation (1), which is line 2 in the equation system.