Solving coupled set of PDEs and ODEs for 1D-problem using pde1dM
17 views (last 30 days)
Show older comments
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
0 Comments
Answers (1)
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).
See Also
Categories
Find more on Geometry and Mesh in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!