필터 지우기
필터 지우기

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 CenterFile Exchange에서 Equation Solving에 대해 자세히 알아보기

제품


릴리스

R2021a

Community Treasure Hunt

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

Start Hunting!

Translated by