Index: mttroot/mtt/bin/mtt ================================================================== --- mttroot/mtt/bin/mtt +++ mttroot/mtt/bin/mtt @@ -13,10 +13,13 @@ ############################################################### ## Version control history ############################################################### ## $Header$ ## $Log$ +## Revision 1.293.2.6 2001/03/03 06:50:38 geraint +## Added #SUMMARY lines for ode2odes. +## ## Revision 1.293.2.5 2001/03/03 00:27:14 geraint ## Fixed ar options to work with GNU ar. Allow mtt to create dependencies for mtt_%.cc when making mtt_%.oct. ## ## Revision 1.293.2.4 2001/03/02 00:45:21 geraint ## Separated Euler and Implicit methods in .cc code and dependencies. @@ -1081,12 +1084,16 @@ ;; implicit) integration_method=implicit; mtt_switches="$mtt_switches implicit"; ;; + rk4) + integration_method=rk4; + mtt_switches="$mtt_switches rk4"; + ;; *) - echo $1 is an unknown integration method - use euler or implicit; + echo $1 is an unknown integration method - use euler, rk4 or implicit; exit;; esac;; -s ) mtt_switches="$mtt_switches $1"; sensitivity=sensitivity ;; @@ -1219,11 +1226,11 @@ echo ' -abg start at abg.m representation' echo ' -c c-code generation' echo ' -d use directory ' echo ' -dc Maximise derivative (not integral) causality' echo ' -dc Maximise derivative (not integral) causality' - echo ' -i Use implicit or euler integration' + echo ' -i Use implicit, euler or rk4 integration' echo ' -o ode is same as dae' echo ' -oct use oct files in place of m files where appropriate' echo ' -opt optimise code generation' echo ' -p print environment variables' echo ' -partition partition hierachical system' @@ -2094,16 +2101,16 @@ @echo "Creating $1_ode2odes_common.o" ar -cr \$@ \$^ $1_ode2odes_common.oct : $1_input.oct $1_logic.oct $1_numpar.oct $1_simpar.oct $1_state.oct @echo > /dev/null -$1_ode2odes_euler.m : $1_ode.m $1_odeo.m +$1_ode2odes_euler.m $1_ode2odes_rk4.m : $1_ode.m $1_odeo.m @echo > /dev/null -$1_ode2odes_euler.cc : $1_ode.cc $1_odeo.cc +$1_ode2odes_euler.cc $1_ode2odes_rk4.cc : $1_ode.cc $1_odeo.cc @echo > /dev/null -$1_ode2odes_euler.o : $1_ode.o $1_odeo.o mtt_euler.o - @echo "Creating $1_ode2odes_euler.o" +$1_ode2odes_euler.o $1_ode2odes_rk4.o : $1_ode.o $1_odeo.o mtt_euler.o + @echo "Creating \$@" ar -cr \$@ \$^ $1_ode2odes_euler.oct : $1_ode.oct $1_odeo.oct mtt_euler.oct @echo > /dev/null $1_ode2odes_implicit.m : $1_cseo.m $1_csex.m $1_smxa.m $1_smxax.m @@ -2413,11 +2420,11 @@ mtt $mtt_switches -q -u $1 ode2odes oct else make_ode2odes $1 m $integration_method endif endif -ifeq ($integration_method,euler) +ifneq ($integration_method,implicit) $1_ode2odes.m : $1_def.r $1_sympars.txt\ $1_simpar.m $1_numpar.m $1_state.m $1_input.m \ $1_ode.m $1_odeo.m $1_logic.m ifeq ($using_oct,yes) echo "*** Warning: Shouldn't be here! Creating dummy $1_ode2odes.m" @@ -2496,11 +2503,11 @@ $1_ode2odes.p : $1_ode2odes.m $1_def.r $1_smxa.p $1_smxax.p\ $1_simpar.p $1_numpar.p $1_state.p $1_input.p \ $1_csex.p $1_cseo.p $1_logic.p mtt_m2p $1_ode2odes.m $integration_method $stdin endif -ifeq ($integration_method,euler) +ifneq ($integration_method,implicit) $1_ode2odes.p : $1_ode2odes.m $1_def.r\ $1_simpar.p $1_numpar.p $1_state.p $1_input.p \ $1_ode.p $1_odeo.p $1_logic.p mtt_m2p $1_ode2odes.m $integration_method $stdin endif Index: mttroot/mtt/bin/trans/make_ode2odes ================================================================== --- mttroot/mtt/bin/trans/make_ode2odes +++ mttroot/mtt/bin/trans/make_ode2odes @@ -7,10 +7,13 @@ ############################################################### ## Version control history ############################################################### ## $Id$ ## $Log$ +## Revision 1.51.2.6 2001/03/16 03:56:13 geraint +## Removed psignal/siginfo.h - problematic and unnecessary. +## ## Revision 1.51.2.5 2001/03/12 23:16:37 geraint ## Minor improvements to signal handling (.exe). ## ## Revision 1.51.2.4 2001/03/12 03:59:30 geraint ## SIGINT (C-c C-c) now causes simulation data to be dumped to MTT.core. @@ -254,12 +257,25 @@ [u] = ${sys}_input(x,y,t,par); # Input [y] = ${sys}_$odeo(x,u,t,par); # Output if mttj==0 mtt_write(t,x,y,$Nx,$Ny); # Write it out endif +EOF + +if [ "$method" = "rk4" ]; then +cat << EOF >> $filename + [k1] = ddt * ${sys}_${ode}(x,u,t,par); + [k2] = ddt * ${sys}_${ode}(x+k1/2,u,t+ddt/2,par); + [k3] = ddt * ${sys}_${ode}(x+k2/2,u,t+ddt/2,par); + [k4] = ddt * ${sys}_${ode}(x+k3,u,t+ddt,par); + [dx] = [k1 + 2.0 * [k2 + k3] + k4] / (6.0 * ddt); +EOF +else +cat << EOF >> $filename [dx] = ${sys}_$ode(x,u,t,par); # State derivative EOF +fi if [ "$method" = "implicit" ]; then cat<< EOF >> $filename [AA] = ${sys}_smxa(x,u,ddt,par); # (I-Adt) and (I-Adt)x @@ -782,12 +798,41 @@ u = mtt_input (x, y, t, par); if (0 == j) { mtt_write (t, x, y, nrows); } +EOF +if [ "$method" = "rk4" ]; then +cat << EOF >> $filename + { + static ColumnVector + k1 (MTTNX,0.0), + k2 (MTTNX,0.0), + k3 (MTTNX,0.0), + k4 (MTTNX,0.0); + + const double + t1 = t + ddt/2.0, + 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 <> $filename AA = mtt_smxa (x, u, ddt, par);