現在この質問をフォロー中です
- フォローしているコンテンツ フィードに更新が表示されます。
- コミュニケーション基本設定に応じて電子メールを受け取ることができます。
How to numerically solve a system of coupled partial differential and algebraic equations?
14 ビュー (過去 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:
Solve for C1, C2, C3, φ.
χ1 -χ7 are known constants.
Initial Conditions:
Boundary conditions
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!エラーが発生しました
ページに変更が加えられたため、アクションを完了できません。ページを再度読み込み、更新された状態を確認してください。
Web サイトの選択
Web サイトを選択すると、翻訳されたコンテンツにアクセスし、地域のイベントやサービスを確認できます。現在の位置情報に基づき、次のサイトの選択を推奨します:
また、以下のリストから Web サイトを選択することもできます。
最適なサイトパフォーマンスの取得方法
中国のサイト (中国語または英語) を選択することで、最適なサイトパフォーマンスが得られます。その他の国の 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 (한국어)