フィルターのクリア

Why do I receive an error in calling "decic"?

6 ビュー (過去 30 日間)
Tobias Hubach
Tobias Hubach 2022 年 1 月 26 日
回答済み: Ishu 2023 年 11 月 28 日
Hey everyone,
I'am trying to solve a DAE-system with the help of ode15s. I already tried ode15i in another code, but had some problems. So I want to check if ode15s performs better.
When computing the consistent initial conditions for ode15i I receive the following error message:
Not enough input arguments.
Error in DGL_Membrane>@(t,y,yp)M(t,y)*yp-F(t,y) (line 101)
implicitDAE = @(t,y,yp) M(t,y)*yp - F(t,y)
Error in decic (line 66)
res = feval(odefun,t0,y0,yp0,varargin{:});
Error in DGL_Membrane (line 102)
[y0,yp0] = decic(implicitDAE,0,y0est,[],yp0est,[],opt);
My code is the following (following the procedure of https://de.mathworks.com/help/symbolic/solve-daes-using-mass-matrix-solvers.html)
function [ySol,tSol] = DGL_Membrane(a,k,L)
syms asym bsym csym dsym fsym hsym isym ksym jsym x(t) y(t) z(t) xx(t) e(t)
Dx = diff(x);
Dy = diff(y);
Dz = diff(z);
Dxx = diff(xx);
Eqn1 = -bsym * diff(x(t),1) - x*bsym*csym*dsym*diff(y(t),1) == asym;
Eqn2 = -fsym*diff(z(t),1) - z * fsym * hsym * dsym * diff(y(t),1) == e;
Eqn3 = -isym*Dxx(t) - xx * isym*jsym * dsym * diff(y(t),1) == ksym;
Eqn4 = csym * x(t) + hsym* z(t) + jsym*xx(t) ==0;
Eqn5 = csym * asym + hsym* e(t) + jsym*ksym == 0;
DAEs=[Eqn1,Eqn2,Eqn3,Eqn4,Eqn5];
DAEvars = [x(t);y(t);z(t);xx(t);e(t)];
[M,f] = massMatrixForm(DAEs,DAEvars);
pDAEs = symvar(DAEs);
pDAEvars = symvar(DAEvars);
extraParams = setdiff(pDAEs, pDAEvars);
M = odeFunction(M, DAEvars,bsym,csym,dsym,fsym,hsym,isym,jsym);
f = odeFunction(f, DAEvars,asym,csym,hsym,jsym,ksym);
% Declare parameters
asym = a;
bsym=1e-6*L; %P1
csym=-1; %Z1
dsym=9.64853321233100184e4/(8.3145*293.15); %F/RT
fsym=0.01e-6*L; %P2
hsym=-2;%Z2
isym=10e-6*L; %P3
jsym=1; %Z3
ksym=k;
c1f=0.2e-3; %feed concentrations mol/m^3
c2f=0.2e-3;
c3f=0.6e-3;
F = @(t,Y) f(t,Y,asym,csym,hsym,jsym,ksym);
y0est = [c1f; 0; c2f;c3f;(-csym * asym - jsym*ksym)/hsym];
yp0est = zeros(5,1);
opt = odeset('Mass', M, 'InitialSlope', yp0est,...
'RelTol', 10.0^(-6), 'AbsTol' , 10.0^(-6));
implicitDAE = @(t,y,yp) M(t,y)*yp - F(t,y)
[y0,yp0] = decic(implicitDAE,0,y0est,[],yp0est,[],opt);
opt = odeset(opt, 'InitialSlope', yp0);
time=linspace(0,L,500);
[tSol,ySol] = ode15s(F,time,y0est,opt);
plot(tSol,ySol(:,[1,3,4,5]),'LineWidth',2)
end
Does anyone have an idea what's wrong with my decic-input?
Thanks for help and best regards !
  3 件のコメント
Tobias Hubach
Tobias Hubach 2022 年 1 月 26 日
Hello Torsten,
When I refer to the decic documentation it says 7 inputs (or 6 without options)
Torsten
Torsten 2022 年 1 月 26 日
編集済み: Torsten 2022 年 1 月 26 日
You are right, it was the symbolic decic function.
I think you will have to define MM the same as F
MM = @(t,y) M(t,y,asym,csym,hsym,jsym,ksym)
and use MM in the definition of implicitDAE.

サインインしてコメントする。

回答 (1 件)

Ishu
Ishu 2023 年 11 月 28 日
Hi Tobias Hubach,
I understand that you are trying to use the "decic" function of MATLAB to solve the differential equations.
I tried to reproduce the error at my end and it also shows me the same error. After that I checked the function handles "M" and "F" and the input arguments required by the handle "M" are a lot not just the "(t,y)". Because of this reason it gives the error of "Not enough input arguments".
Further I looked into the code and found that you are changing the defination of "M" after its initialization. Explained in the given code snippet.
[M,f] = massMatrixForm(DAEs,DAEvars); % you have initialised here
M
M = 
pDAEs = symvar(DAEs);
pDAEvars = symvar(DAEvars);
extraParams = setdiff(pDAEs, pDAEvars);
M = odeFunction(M, DAEvars,bsym,csym,dsym,fsym,hsym,isym,jsym) % here you have changed "M"
M = function_handle with value:
@(t,in2,param1,param2,param3,param4,param5,param6,param7)reshape([-param1,0.0,0.0,0.0,0.0,-param1.*param2.*param3.*in2(1,:),-param3.*param4.*param5.*in2(3,:),-param3.*param6.*param7.*in2(4,:),0.0,0.0,0.0,-param4,0.0,0.0,0.0,0.0,0.0,-param6,0.0,0.0,0.0,0.0,0.0,0.0,0.0],[5,5])
Hope it helps.

カテゴリ

Help Center および File ExchangeEquation Solving についてさらに検索

製品


リリース

R2021a

Community Treasure Hunt

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

Start Hunting!

Translated by