How to numerically solve a system of coupled partial differential and algebraic equations?
19 次查看(过去 30 天)
显示 更早的评论
Arpita
2023-1-25
I have a system of coupled partial differential and algebraic equations.
Two 1-D parabolic pdes coupled (function of x and time) with two algebraic equations. What would be the way to solve this sytem?
14 个评论
Torsten
2023-1-25
Without knowing your equations together with initial and boundary conditions, nothing useful can be said.
Arpita
2023-1-26
Hi Torsten,
I really appreciate you getting back to me.
There is a possibility of singularity in the matrix. My equations include:
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1275045/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1275050/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1275145/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1275060/image.png)
Solve for C1, C2, C3, φ.
χ1 -χ7 are known constants.
Initial Conditions:
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1275065/image.png)
Boundary conditions
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1275070/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1275075/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1275080/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1275085/image.png)
Torsten
2023-1-26
编辑:Torsten
2023-1-26
This will be hard to solve, and there is no MATLAB code that could help you.
The usual procedure is to discretize the spatial derivatives in equations (1) and (2) and solve the resulting system of differential-algebraic equations using ODE15S.
But you'll have a lot of trouble with it, that's for sure.
Arpita
2023-1-26
Would that something similar to what's done in this example:
https://www.mathworks.com/matlabcentral/answers/1730075-how-to-solve-coupled-partial-differential-equations-with-method-of-lines
Torsten
2023-1-26
编辑:Torsten
2023-1-26
Yes, the method-of-lines means to use difference quotients as approximations for the spatial derivatives (d/dx, d^2/dx^2) and transform the partial differential-algebraic system into an ordinary differential-algebraic system.
This can then be solved using a MATLAB solver like ODE15S.
Arpita
2023-1-26
I convert my pdes into algebraic equations (finite difference method) as also done in the example. Do I just add the algebraic equations into the for loop?
Torsten
2023-1-26
Yes, but you must tell the solver that these are algebraic and not differential equations.
This is done by zero rows within the mass matrix M which is part of the differential-algebraic system in the form of
M*y' = f(t,y)
Arpita
2023-2-2
In solving these coupled PDEs and algebraic equations I am facing the following error:
"Warning: Failure at t=1.048722e-05. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (3.725809e-20) at time t."
I tried changing the solver to ode23t to skirt around stiffness issues, if any, but that doesn't solve the problem.
Setting 'Mass Singular" to 'yes" does not help either.
Could someone help with this problem?
Torsten
2023-2-2
This means that ODE15S doesn't like your DAE system.
The reasons can be manifold: programming errors, difficulty of the system to be integrated, ...
So no advice can be given here.
Arpita
2023-2-2
编辑:Torsten
2023-2-2
I don't think I have any programming errors, but could you possibly help me?
If the issue is the difficulty of the system to be integrated, should I move to a different software?
clear
clc
%% Parameters %%
% Experimental conditions %
p.Q1=37.33/60;% 4.4664/60; % flowrate to the stack (mL/s)
p.I=-1258/1120*10^(-3); % applied current density (A/cm2)
p.T1=500;% round(41/Q1); % charging time(s)
p.C0=0.0342; % initial concentration (mmol/ml)
% Electrode and cell physical properties %
p.Mass=10; %electrode weight (g)
p.N=1; % number of electrode pairs
p.psp=0.708; % spacer porosity
p.pma=0.4; % macroprosity
p.pmi=0.28; % microporosity
p.Lelec=0.028; % each electrode thickness (cm)
p.Lsp=0.01; % spacer thickness (cm)
p.a=1120; % each electrode area (cm^2)
p.hrt1=p.Lsp*p.a*p.N/p.Q1; % hydraulic retention time (s) in the spacer
p.tlag=round(1*p.hrt1);
p.alpha=50;
p.Rext=200;% resistance (om*cm2)
% Electrochemical properties %
p.De=1.68*10^(-5); % effective diffusion coefficient (cm^2/s)
p.R=6.1;% specific electrode resistance (om*mmol/cm), from zhao, WR, 'optimation' 2013
p.u=0; % chemical attraction term (kT)
p.Vt=0.0256; % thermal voltage (V)
p.Cst1=72; % stern capacitance (F/mL) value from JE.Dykstra W 2016.
p.F=96.5; % farady constant (C/mmol)
p.i=p.I/p.F; % convert (A/cm2) to (mmol/s/cm2)
p.Vapplied = (1.2/2/p.Vt);
p.V1 = p.Vapplied-(p.i*p.Lsp/(4*p.C0*p.psp*p.De));
% using method of lines %
tf = 200; % End time
p.nz = 100; % Number of nodes
x = linspace(0,p.Lelec,p.nz); % Space domain
dx = diff(x);
p.dx = dx(1); % Space increment
%% Initial conditions %%
Ceff_0 = p.C0*(p.pma+p.pmi); % Ceff(x,0)= C0
Cch_0 = 0; % Cch(x,0) = 0
Cma_0 = p.C0; % Cma(x,0) = C0
phima_0 = p.V1; %phima(x,0) = V1
Csp = p.C0.*ones(p.nz,1);
Ceff = ones(p.nz,1);
Ceff=Ceff_0.*Ceff;
Cch = zeros(p.nz,1);
Cma = ones(p.nz,1);
Cma=Cma_0.*Cma;
phima = ones(p.nz,1);
phima=phima_0.*phima;
%y0 = [Csp;Ceff;Cch;Cma;phima];
y0 = [Ceff;Cch;Cma;phima];
tspan = [0 tf];
%% Boundary conditions %%
dCmadx = 0; % dCmadx(Lelec,t) = 0
dphimadx =0; % dphimadx(Lelec,t) =0
%% Mass matrix (for DAE system) %%
M = eye(4*p.nz);
for i=(2*p.nz+1):(4*p.nz)
M(i,i)=0;
end
options = odeset('Mass',M,'RelTol',1e-3,'AbsTol',1e-10);%'MassSingular','yes');
[t,y] = ode23t(@(t,y)porada1D(t,y,p),tspan,y0,options);
%% Plot %%
Ceff_matrix=y(:,1:p.nz);
Cch_matrix=y(:,p.nz+1:2*p.nz);
Cma_matrix=y(:,2*p.nz+1:3*p.nz);
phima_matrix=y(:,3*p.nz+1:4*p.nz);
plot(t,Ceff_matrix(:,:))
grid
xlabel('time')
ylabel('Ceff')
%title('Dimensionless concentration at ''z = L''')
%figure()
%plot(t, phima_matrix(:,:))
%% function %%
function out = porada1D(~,y,p)
% Variables
Ceff = y(1:p.nz,1);
Cch = y(p.nz+1:(2*p.nz),1);
Cma = y((2*p.nz)+1:(3*p.nz),1);
phima = y((3*p.nz)+1:(4*p.nz),1);
%phima = y((4*p.nz)+1:(5*p.nz),1);
% Space derivatives
dCmadx = zeros(p.nz,1);
d2Cmadx2 = zeros(p.nz,1);
dphimadx = zeros(p.nz,1);
d2phimadx2 = zeros(p.nz,1);
dCmadx(2:end-1) = (Cma(3:end)-Cma(1:end-2))./(2.*p.dx);
d2Cmadx2(2:end-1) = (Cma(3:end)-(2*Cma(2:end-1))+Cma(1:end-2))./(p.dx^2);
dphimadx(2:end-1) = (phima(3:end)-phima(1:end-2))./(2.*p.dx);
d2phimadx2(2:end-1) = (phima(3:end)-(2*phima(2:end-1))+Cma(1:end-2))./(p.dx^2);
% Governing equations
%dCspdt = (2.*p.pma.*p.De/p.psp/p.Lsp).*dCmadx + ((1./p.psp./p.hrt1).*(p.C0-Csp));
dCeffdt = p.pma.*p.De.*d2Cmadx2;
dCchdt = 2.*(p.pma./p.pmi).*p.De.*(dCmadx.*dphimadx+Cma.*d2phimadx2);
phit = phima+(asinh(Cch./(-2.*Cma)))-(p.F.*Cch./(p.Vt.*(p.Cst1+p.alpha.*(Cch.*Cch))))-p.V1;
Cefft = Ceff-p.pma.*Cma-(p.pmi.*Cma.*cosh(Cch./(-2.*Cma)));
% Boundary conditions
Cma(1,1)=p.C0;
dCmadx(end)=0;
dphimadx(end)=0;
%out = [dCspdt;dCeffdt;dCchdt;phit;Cefft];
out = [dCeffdt;dCchdt;phit;Cefft];
end
Torsten
2023-2-2
% Boundary conditions
Cma(1,1)=p.C0;
dCmadx(end)=0;
dphimadx(end)=0;
The boundary conditions you set at the end of your function "porada1D" are nowhere used in the calculation of your output vector
out = [dCeffdt;dCchdt;phit;Cefft];
This cannot be correct.
Arpita
2023-2-3
The version I sent you had the boundary conditions after the governing equations. Sorry about that.
Even with boundary conditions before the governing equations, I see the same problem.
With all three boundary conditions, I get the error that the DAE is greater than 1.
When I only use the dCmadx = 0 and dphimadx = 0 as the boundary conditions, I get the error: " Unable to meet integration tolerances without reducing the step size below the smallest value allowed (4.336809e-19) at time t."
Torsten
2023-2-3
You should start with a simpler problem from which you know that potentially arising problems with the integrator stem from your programming, not from the difficulty or even unsolvability of the problem itself.
回答(1 个)
Sarthak
2023-3-9
Hi,
One way to solve a system of coupled partial differential equations (PDEs) and algebraic equations is to use a numerical method such as finite difference or finite element method. Here is an outline of the steps involved:
- Discretize the system of PDEs using a numerical method such as finite difference or finite element method. This will transform the PDEs into a system of algebraic equations.
- Combine the discretized PDEs with the algebraic equations to form a system of nonlinear algebraic equations.
- Use a numerical solver such as the Newton-Raphson method or a quasi-Newton method to solve the system of nonlinear algebraic equations.
- Repeat the process for each time step to obtain a time-dependent solution.
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Boundary Conditions 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!发生错误
由于页面发生更改,无法完成操作。请重新加载页面以查看其更新后的状态。
您也可以从以下列表中选择网站:
如何获得最佳网站性能
选择中国网站(中文或英文)以获得最佳网站性能。其他 MathWorks 国家/地区网站并未针对您所在位置的访问进行优化。
美洲
- América Latina (Español)
- Canada (English)
- United States (English)
欧洲
- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)
- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom (English)
亚太
- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)