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.57.2.1 2001/05/04 04:07:24 geraint
## Numerical solution of algebraic equations.
## sys_ae.cc written for unsolved inputs.
## Solution of equations using hybrd from MINPACK (as used by Octave fsolve).
##
## Revision 1.57 2001/04/01 03:38:54 geraint
## Reset row to zero after write to file, ready for subsequent runs.
|
>
>
>
|
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.57.2.2 2001/06/05 03:20:40 geraint
## added -ae option to select algebraic equation solution method.
##
## Revision 1.57.2.1 2001/05/04 04:07:24 geraint
## Numerical solution of algebraic equations.
## sys_ae.cc written for unsolved inputs.
## Solution of equations using hybrd from MINPACK (as used by Octave fsolve).
##
## Revision 1.57 2001/04/01 03:38:54 geraint
## Reset row to zero after write to file, ready for subsequent runs.
|
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
|
args (0) = octave_value (x);
f = feval ("${sys}_state", args, 1);
return f(0).${vector_value} ();
#endif
}
inline ColumnVector
mtt_${ode} (ColumnVector &x,
ColumnVector &u,
const double &t,
ColumnVector &par)
{
#ifdef STANDALONE
return F${sys}_${ode} (x, u, t, par);
#else
static octave_value_list args, f;
args (0) = octave_value (x);
args (1) = octave_value (u);
args (2) = octave_value (t);
args (3) = octave_value (par);
f = feval ("${sys}_${ode}", args, 1);
return f(0).${vector_value} ();
#endif
}
inline ColumnVector
mtt_${odeo} (ColumnVector &x,
ColumnVector &u,
const double &t,
ColumnVector &par)
{
#ifdef STANDALONE
return F${sys}_${odeo} (x, u, t, par);
#else
|
|
|
|
|
|
|
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
|
args (0) = octave_value (x);
f = feval ("${sys}_state", args, 1);
return f(0).${vector_value} ();
#endif
}
inline ColumnVector
mtt_rate (ColumnVector &x,
ColumnVector &u,
const double &t,
ColumnVector &par)
{
#ifdef STANDALONE
return F${sys}_${ode} (x, u, t, par);
#else
static octave_value_list args, f;
args (0) = octave_value (x);
args (1) = octave_value (u);
args (2) = octave_value (t);
args (3) = octave_value (par);
f = feval ("${sys}_${ode}", args, 1);
return f(0).${vector_value} ();
#endif
}
inline ColumnVector
mtt_output (ColumnVector &x,
ColumnVector &u,
const double &t,
ColumnVector &par)
{
#ifdef STANDALONE
return F${sys}_${odeo} (x, u, t, par);
#else
|
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
|
{
x (i) = x0 (i);
}
for (register int j = 0, i = 1; i <= ilast; i++)
{
u = mtt_input (x, y, t, par);
y = mtt_${odeo} (x, u, t, par);
if (0 == j)
{
mtt_write (t, x, y, nrows);
}
EOF
if [ "$method" = "rk4" ]; then
cat << EOF >> $filename
|
|
|
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
|
{
x (i) = x0 (i);
}
for (register int j = 0, i = 1; i <= ilast; i++)
{
u = mtt_input (x, y, t, par);
y = mtt_output (x, u, t, par);
if (0 == j)
{
mtt_write (t, x, y, nrows);
}
EOF
if [ "$method" = "rk4" ]; then
cat << EOF >> $filename
|
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
|
t2 = t + ddt;
ColumnVector
x1 (x),
x2 (x),
x3 (x);
k1 = ddt * mtt_${ode} (x , u, t , par); x1 += k1 * 0.5;
k2 = ddt * mtt_${ode} (x1, u, t1, par); x2 += k2 * 0.5;
k3 = ddt * mtt_${ode} (x2, u, t1, par); x3 += k3;
k4 = ddt * mtt_${ode} (x3, u, t2, par);
dx = (k1 + 2.0 * (k2 + k3) + k4) / (6.0 * ddt);
}
EOF
else
cat << EOF >> $filename
dx = mtt_${ode} (x, u, t, par);
EOF
fi
if [ "$method" = "implicit" ]; then
cat <<EOF >> $filename
AA = mtt_smxa (x, u, ddt, par);
|
|
|
|
|
|
|
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
|
t2 = t + ddt;
ColumnVector
x1 (x),
x2 (x),
x3 (x);
k1 = ddt * mtt_rate (x , u, t , par); x1 += k1 * 0.5;
k2 = ddt * mtt_rate (x1, u, t1, par); x2 += k2 * 0.5;
k3 = ddt * mtt_rate (x2, u, t1, par); x3 += k3;
k4 = ddt * mtt_rate (x3, u, t2, par);
dx = (k1 + 2.0 * (k2 + k3) + k4) / (6.0 * ddt);
}
EOF
else
cat << EOF >> $filename
dx = mtt_rate (x, u, t, par);
EOF
fi
if [ "$method" = "implicit" ]; then
cat <<EOF >> $filename
AA = mtt_smxa (x, u, ddt, par);
|