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

17 views (last 30 days)
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

Answers (1)

Bill Greene
Bill Greene on 1 Feb 2022
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).
  1 Comment
Paul Maurerlehner
Paul Maurerlehner on 1 Feb 2022
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.

Sign in to comment.

Products


Release

R2020a

Community Treasure Hunt

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

Start Hunting!