1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
|
#! /bin/sh
######################################
##### Model Transformation Tools #####
######################################
###############################################################
## Version control history
###############################################################
## $Id$
## $Log$
## Revision 1.63 2001/08/07 04:39:24 geraint
## Consolidated dassl and residual functions.
##
## Revision 1.62 2001/08/01 22:14:32 geraint
## Bug fix for dassl.
##
## Revision 1.61 2001/08/01 04:06:07 geraint
|
>
>
>
|
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
|
#! /bin/sh
######################################
##### Model Transformation Tools #####
######################################
###############################################################
## Version control history
###############################################################
## $Id$
## $Log$
## Revision 1.64 2001/08/08 02:15:00 geraint
## Rationalisation of solver code, beginning with algebraic solvers.
##
## Revision 1.63 2001/08/07 04:39:24 geraint
## Consolidated dassl and residual functions.
##
## Revision 1.62 2001/08/01 22:14:32 geraint
## Bug fix for dassl.
##
## Revision 1.61 2001/08/01 04:06:07 geraint
|
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
|
#endif
}
#ifdef STANDALONE
ColumnVector
Fmtt_residual (const ColumnVector &X,
const ColumnVector &DX,
double t)
{
#else // !STANDALONE
DEFUN_DLD (mtt_residual, args, , "mtt_residual")
{
static ColumnVector X (MTTNX+MTTNYZ);
static ColumnVector DX (MTTNX+MTTNYZ);
static double t;
X = args(0).${vector_value} ();
DX = args(1).${vector_value} ();
t = args(2).double_value ();
#endif // STANDALONE
static ColumnVector residual (MTTNX+MTTNYZ);
static ColumnVector U (MTTNU+MTTNYZ);
static ColumnVector u (MTTNU);
static ColumnVector y (MTTNY,0.0);
static ColumnVector par (MTTNPAR);
|
|
>
>
>
|
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
|
#endif
}
#ifdef STANDALONE
ColumnVector
Fmtt_residual (const ColumnVector &X,
const ColumnVector &DX,
double t,
int &ires)
{
#else // !STANDALONE
DEFUN_DLD (mtt_residual, args, , "mtt_residual")
{
static ColumnVector X (MTTNX+MTTNYZ);
static ColumnVector DX (MTTNX+MTTNYZ);
static double t;
static int &ires;
X = args(0).${vector_value} ();
DX = args(1).${vector_value} ();
t = args(2).double_value ();
ires = static_cast<int>(args(3).double_value ());
#endif // STANDALONE
static ColumnVector residual (MTTNX+MTTNYZ);
static ColumnVector U (MTTNU+MTTNYZ);
static ColumnVector u (MTTNU);
static ColumnVector y (MTTNY,0.0);
static ColumnVector par (MTTNPAR);
|
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
|
yz = feval ("${sys}_ae", new_args, 1)(0).${vector_value} ();
#endif
for (register int i = 0; i < MTTNX; i++)
residual (i) = dx(i) - DX(i);
if (MTTNYZ > 0)
residual.insert (yz,MTTNX);
|
>
|
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
|
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
|
yz = feval ("${sys}_ae", new_args, 1)(0).${vector_value} ();
#endif
for (register int i = 0; i < MTTNX; i++)
residual (i) = dx(i) - DX(i);
if (MTTNYZ > 0)
{
residual.insert (yz,MTTNX);
// XXX:
// DASSL needs residual to be dependent on Ui and Uidot
// mtt_dassl always sets the initial Ui to zero, so
// Ui - h*Uidot should equal zero BUT, we don't know h
static double t_old;
double step;
if (t != t_old)
{
step = t - t_old;
t = t_old;
}
else
step = t;
for (register int i = MTTNX; i < MTTNX+MTTNYZ; i++)
residual(i) += X(i) - DX(i)*step;
}
#ifdef STANDALONE
return residual;
#else // !STANDALONE
return octave_value (residual);
#endif // STANDALONE
}
|