ODE15s / DAE error message
Daniel
I'm trying to solve a DAE system of the following form:
M=diag([1 1 1 1 1 1 0 0 0 0]);
options=odeset('Mass',M);
[t,sol]=ode15s(@fun,tspan,x0,options)
function d_dt=fun(t,x)
c1=... some expression
c2=...
c3=...
c4=...
c5=...
c6=...
res7=c7+c8
res8=c1*c7
res9=parameter/(c1+c2+c3)*c7
res10=parameter/(c1+c2+c3)*c8
d_dt=[c1 c2 c3 c4 c5 c6 res7 res8 res9 res10 ]
I get the error message:
###################################################################
??? Error using ==> mtimes
Inner matrix dimensions must agree.
Error in ==> daeic12 at 62
F = UM' * f;
Error in ==> ode15s at 395
[y,yp,f0,dfdy,nFE,nPD,Jfac] = daeic12(odeFcn,odeArgs,t,ICtype,Mt,y,yp0,f0,...
###################################################################
So question 1:
What is this error about???
question 2:
What does the option 'Vectorized', 'On' mean and could it help in this case?
Thanks for your help in advance!
Steven_Lord
"Daniel " <ddr_l...@gmx.net> wrote in message
news:i8co5t$192$1...@fred.mathworks.com...
> Hi there,
>
> I'm trying to solve a DAE system of the following form:
>
> M=diag([1 1 1 1 1 1 0 0 0 0]);
> options=odeset('Mass',M);
> [t,sol]=ode15s(@fun,tspan,x0,options)
>
> function d_dt=fun(t,x)
> c1=... some expression
> c2=...
> c3=...
> c4=...
> c5=...
> c6=...
> res7=c7+c8
> res8=c1*c7
> res9=parameter/(c1+c2+c3)*c7
> res10=parameter/(c1+c2+c3)*c8
>
> d_dt=[c1 c2 c3 c4 c5 c6 res7 res8 res9 res10 ]
The ODE solvers expect your function to return a _column_ vector of values,
not a _row_ vector.
http://www.mathworks.com/help/techdoc/ref/ode23.html
"Function f = odefun(t,y), for a scalar t and a column vector y, must return
a column vector f corresponding to f(t, y)."
*snip*
> So question 1: What is this error about???
Well, you haven't included your full code, but I suspect the different
orientations is the cause of this particular error.
> question 2:
> What does the option 'Vectorized', 'On' mean and could it help in this
> case?
The various options available for use with the ODE solvers are described
here:
http://www.mathworks.com/help/techdoc/ref/odeset.html#f92-1023074
I don't think using that option would necessarily resolve that particular
error.
--
Steve Lord
sl...@mathworks.com
comp.soft-sys.matlab (CSSM) FAQ: http://matlabwiki.mathworks.com/MATLAB_FAQ
To contact Technical Support use the Contact Us link on
http://www.mathworks.com
Daniel
> > d_dt=[c1 c2 c3 c4 c5 c6 res7 res8 res9 res10 ]
>
> The ODE solvers expect your function to return a _column_ vector of values,
> not a _row_ vector.
>
> http://www.mathworks.com/help/techdoc/ref/ode23.html
>
> "Function f = odefun(t,y), for a scalar t and a column vector y, must return
> a column vector f corresponding to f(t, y)."
>
Sorry, just a typing error. The d_dt vector is a column vector, I checked that.
So any further ideas what might cause this error?
(Shall I give the full code?)
Torsten Hennig
>
> I'm trying to solve a DAE system of the following
> form:
>
> M=diag([1 1 1 1 1 1 0 0 0 0]);
> options=odeset('Mass',M);
> [t,sol]=ode15s(@fun,tspan,x0,options)
>
Did you initialize x0 ?
> function d_dt=fun(t,x)
> c1=... some expression
> c2=...
> c3=...
> c4=...
> c5=...
> c6=...
> res7=c7+c8
> res8=c1*c7
> res9=parameter/(c1+c2+c3)*c7
> res10=parameter/(c1+c2+c3)*c8
>
> d_dt=[c1 c2 c3 c4 c5 c6 res7 res8 res9 res10 ]
>
> I get the error message:
>
> ######################################################
> #############
> ??? Error using ==> mtimes
> Inner matrix dimensions must agree.
>
> Error in ==> daeic12 at 62
> F = UM' * f;
>
> Error in ==> ode15s at 395
> [y,yp,f0,dfdy,nFE,nPD,Jfac] =
> c] = daeic12(odeFcn,odeArgs,t,ICtype,Mt,y,yp0,f0,...
> ######################################################
> #############
>
> So question 1:
> What is this error about???
>
> question 2:
> What does the option 'Vectorized', 'On' mean and
> could it help in this case?
>
> Thanks for your help in advance!
Best wishes
Torsten.
Daniel
> Did you initialize x0 ?
Yes I did!
As I understand from the help, 'Vectorize', 'On' is applicable for this problem.
Unfortunately this results in another (quiet similar) error message:
#####################################################################
??? Error using ==> minus
Matrix dimensions must agree.
Error in ==> odenumjac at 148
Fdiff = Fdel - Fvalue(:,ones(1,ny));
Error in ==> daeic12 at 38
[DfDy,Joptions.fac,nF] = odenumjac(fun, {t0,y,args{:}}, f, Joptions);
Error in ==> ode15s at 395
[y,yp,f0,dfdy,nFE,nPD,Jfac] = daeic12(odeFcn,odeArgs,t,ICtype,Mt,y,yp0,f0,...
###################################################################
Daniel
The description for daeic12 says:
' DAEIC12 attempts to find a set of consistent initial conditions for
equations of the form M(t)*y' = f(t,y) '
So what does this mean exactly?
I guess it has something to do with the initials I set...but why consistent? I only can take the initials which are physical meaningful ...
Torsten Hennig
You have four algebraic equations in your DAE-system
(namely the last 4).
If you don't supply starting values for x1,...,x10
such that f7(t=0,x1,...,x10), f8(t=0,x1,...,x10),
f9(t=0,x1,...,x10), f10(t=0,x1,...,x10) are identically
zero, DAEIC12 attempts to adjust your x7,...,x10 such
that this is the case.
But this has nothing to do with the error message you
receive, I would guess.
Best wishes
Torsten.
Daniel
But the iteration does not start because of this error. The value of the t variable stays zero allthough the function is called about 100 times. Then the ode returns this error message caused by the inconsistency in the DAE initial conditions...as far as I understand this. But why do I need to choose the 'right' initial guesses?
Torsten Hennig
The error messages you posted so far said something
about 'matrix dimensions that do not agree', but nothing
about 'an inconsistency in the DAE initial conditions'.
> But why do I need to choose the
> 'right' initial guesses?
If you have an algebraic equation like y(7) - 500 = 0
in your system, you want y(7) = 500 for all times t.
If at time t = 0, you provide y(7) = 2000,
you provide a _wrong_ initial value.
The solver tries to adjust it to y(7) = 500,
but if your initial guess is too far away form the correct value, the solver may fail.
This is different from differential variables
where there is no 'correct' initial condition, but you
can start from where you like.
Best wishes
Torsten.
Daniel
> > Torsten Hennig <Torsten...@umsicht.fhg.de> wrote
> > in message
> > <652704926.53626.1286282352510.JavaMail.root@gallium.m
> > athforum.org>...
> in your system, you want y(7) = 500 for all times t.
> If at time t = 0, you provide y(7) = 2000,
> you provide a _wrong_ initial value.
> The solver tries to adjust it to y(7) = 500,
> but if your initial guess is too far away form the correct value, the solver may fail.
>
> This is different from differential variables
> where there is no 'correct' initial condition, but you
> can start from where you like.
>
> Best wishes
> Torsten.
OK. Thank you Torsten!!!
But what if y (as you introduced it above) is a vector so that I get something like y(1:5)-500=0 ? Is this a case for 'Vectorized', 'On' now?
I arrive at this case for my DAE system from a discretization with the finite volume method. Its kind of difficult to get the initial values for this then...because I know y(1) and maybe y(5), but not the values in between.
Daniel
edit: y(1:5)-500*=0
* - can also be a vector e.g. z(1:5)
Daniel
Found the error which was causing the error messages above: the Mass-Matrix must a Matrix. Stupid mistake but very simple to solve...
But another error occurred now:
'??? Error using ==> daeic12 at 167 Need a better guess y0 for consistent initial conditions.'
Again the question: what shall I do if I do not know all initial conditions? Is there a way to tell the ode solver? Or is there a way to set the conditions physically senseless but at least consistent?
Torsten Hennig
> Found the error which was causing the error messages
> above: the Mass-Matrix must a Matrix. Stupid mistake
> but very simple to solve...
>
> But another error occurred now:
> '??? Error using ==> daeic12 at 167 Need a better
> guess y0 for consistent initial conditions.'
>
> Again the question: what shall I do if I do not know
> all initial conditions? Is there a way to tell the
> ode solver? Or is there a way to set the conditions
> physically senseless but at least consistent?
You have to start from a physically senseful solution
because it is this solution which is the basis
for the following time integration.
You have two possibilities to calculate this solution:
- by hand
- by using MATLAB's fsolve before calling the
DAE-integrator.
Best wishes
Torsten.
Daniel
> > options=odeset('Mass',M);
> > [t,sol]=ode15s(@fun,tspan,x0,options)
> > function d_dt=fun(t,x)
> > c1=... some expression
> > c2=...
> > c3=...
> > c4=...
> > c5=...
> > c6=...
> > res7=c7+c8
> > res8=c1*c7
> > res9=parameter/(c1+c2+c3)*c7
> > res10=parameter/(c1+c2+c3)*c8
> d_dt=[c1 c2 c3 c4 c5 c6 res7 res8 res9 res10 ]
Torsten Hennig <Torsten...@umsicht.fhg.de> wrote in message <559427162.89245.12868...@gallium.mathforum.org>...
OK! Thank you Torsten.
But I have another question! (seems to be a very complicated system that's why that may questions)
I'm referring to the source code above:
For a special region (e.g. x<x1=10m) some of the algebraic terms are not (physically meaningful) present anymore. Is it sufficient to tell Matlab with res8=0 that the equation is not to be considered anymore? I think repeatedly calling the function fun(t,x) with separate values of x is not the solution, isn't it?
Torsten Hennig
> > > options=odeset('Mass',M);
> > > [t,sol]=ode15s(@fun,tspan,x0,options)
>
> > > function d_dt=fun(t,x)
> > > c1=... some expression
> > > c2=...
> > > c3=...
> > > c4=...
> > > c5=...
> > > c6=...
> > > res7=c7+c8
> > > res8=c1*c7
> > > res9=parameter/(c1+c2+c3)*c7
> > > res10=parameter/(c1+c2+c3)*c8
>
> > d_dt=[c1 c2 c3 c4 c5 c6 res7 res8 res9 res10 ]
>
> Torsten Hennig <Torsten...@umsicht.fhg.de> wrote
> in message
The variables corresponding to some algebraic equations
are no longer needed just from the beginning of the
computation or at some stage during the computation ?
> Is it sufficient to tell Matlab with
> res8=0 that the equation is not to be considered
> anymore? I think repeatedly calling the function
> fun(t,x) with separate values of x is not the
> solution, isn't it?
Telling MATLAB res8 = 0 is not a good idea because
the Jacobian matrix will become singular.
Set res8 to a simple algebraic equation that includes
the variable to be determined by the 8th algebraic
equation, e.g. res8 = x(8). Then x(8)=0 for all times t.
Best wishes
Torsten.
Daniel
> > in message
> > <559427162.89245.1286865055825.JavaMail.root@gallium.m
> > athforum.org>...
>
> are no longer needed just from the beginning of the
> computation or at some stage during the computation ?
Again, thank you very much for your help here Torsten!!! I really appreciate it.
Maybe I need to explain the problem a little bit more in specific:
I need to solve a PDE/ODE system for a battery electrode and a battery separator in 1D like shown below:
n= 1 n=6 n=10
|x=0|---------------------------------------------------------------------- |x=x1|-------------------------|x=x2|
< porous electrode > <separator>
solid and solved species only solved species
The porous electrode is filled with solid metal particles and liquid electrolyte solution in between these particles. There are two solid species and 4 solved ones in solution.
Therefore I get 6 species PDEs like
dc/dt =reaction + diffusion + convection + migration = function(c(x), dc/dx, d²c/dx², phi_sol(x)).
The migration is only valid for charged species in solution and is coupled with the 'drag-force' in dependence of the solution potential at the current position x, phi_sol(x).
The potential is calculated with 4 extra first order ODEs in space with the help of the current and potential in solid and solution phase. They are not time dependent.
The DAE system arises because I do a finite volume discretization for the PDE, let's say for n=10 elements...so about n=6 for the electrode and n=4 for the separator. In total there would be a set of 10*10=100 equations.
But I get 6*6=36 ODEs in time for the species and 4*6=24 algebraic equations in the electrode. Consequently, in the separator this is reduced by 2*4=8 ODEs (since the solid components are not present anymore) and 2*4=8 algebraic equations (since the current and potential in the solid metal phase are missing). In total I only need 84 equations. Clear?
So that's why I asked if res8(from x1 to x2)=0 would do?
I also thought of calling separate functions for electrode and separator. But then the time in both functions is different and I would need the boundary values at x=x1 at the exact same time. Right?
Or do you have any other idea how to solve the PDE/ODE system without discretization in Matlab?
Torsten Hennig
> wrote
> > > in message
> > >
> <559427162.89245.1286865055825.JavaMail.root@gallium.m
> > > athforum.org>...
> >
> > The variables corresponding to some algebraic
> equations
> > are no longer needed just from the beginning of the
>
> > computation or at some stage during the computation
> ?
>
> Again, thank you very much for your help here
> Torsten!!! I really appreciate it.
>
>
> Maybe I need to explain the problem a little bit more
> in specific:
>
> I need to solve a PDE/ODE system for a battery
> electrode and a battery separator in 1D like shown
> below:
>
> n= 1
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
> n=6
> n=6
> n=6
> n=6
> n=6
> n=6
> n=6
> n=6
> n=6
> n=6 n=10
> |x=0|-------------------------------------------------
> ---------------------
> |x=x1|-------------------------|x=x2|
> < porous electrode
> rous electrode >
>
>
>
>
>
> <separator>
> solid and solved
> solved species
> only solved
I don't understand why you make life so difficult.
For each node j in the electrode, define a
solution vector
(y_1^(j),...,y_6^(j),y_7^(j),...,y_10^(j))
and for each node j in the seperator, define a
solution vector
(y_1^(j),...,y_6^(j)).
For your above choice of the nuumber of elements,
glue the solution vectors together to an y-vector
of the form
y = (y_1^(1),...,y_10^(1),y_1^(2),...,y_10^(2),
..,y_1^(6),...,y_10^(6),y_1^(7),...,y_6^(7),...,
y_1^(10),...,y_6^(10))
with 60 + 24 = 84 components and call ODE15s
with this vector.
In the function routine, you can write the full
y-vector in local vectors which reflect the underlying
problem (e.g. it is convenient to define subvectors
(y_1^(1),y_1^(2),...,y_1^(10)),
(y_2^(1),y_2^(2),...,y_2^(10)),
..
(y_6^(1),y_6^(2),...,y_6^(10)),
(y_7^(1),y_7^(2),...,y_7^(6)),
..
(y_10^(1),y_10^(2),...,y_10^(6))
for teh different solution variables, I think.)
Best wishes
Torsten.
Daniel
Thanks! Tried that.
But now the DAE index appears to be greater than 1:
??? Error using ==> daeic12 at 77
This DAE appears to be of index greater than 1.
Don't know why exactly because all algebraic states can be described by the algebraic equations...
So again the question can Matlab solve the original PDE/ODE system without discretization?
Torsten Hennig
The reason for this error message usually is a
programming error or inconsistent initial values
for the algebraic variables.
You can check this by calculating the residuals
of the algebraic equations after the first call of
ODE15s to your function routine - they should be very
small if the initial values are consistent.
> So again the question can Matlab solve the original
> PDE/ODE system without discretization?
You may have a look at MATLAB's pdepe, but I think
a mixture of PDEs and ODEs with different domains
on which they are valid requires more flexibility
than pdepe can offer.
Best wishes
Torsten.
dma...@mst.edu
I by chance came across this thread on index error associated with modelling lithium ion battery cell using ode15s in MATLAB. If it will not be too much of a ask, how did you eventually solve this problem because I am presently having this same problem with my battery model using spectral methods.
Regards,
Damola M. Ajiboye.
