| 1 |
/* THIS MATLAB MEX FUNCTION PROVIDES AN INTERFACE TO THE 5TH-ORDER |
|---|
| 2 |
CASH-KARP RUNGE-KUTTA EMBEDDED INTEGRATION SCHEME AS CODED BY |
|---|
| 3 |
NUMERICAL RECIPES IN C, 2ND ED. SEE CHAPTER 16. |
|---|
| 4 |
|
|---|
| 5 |
The Numerical Recipes routines contain the following copyright line: |
|---|
| 6 |
(C) Copr. 1986-92 Numerical Recipes Software 3+)$!|9. |
|---|
| 7 |
|
|---|
| 8 |
The lowest-level routine called is RKCK, which performs one integration |
|---|
| 9 |
over the (adjusted) stepsize. Local error is controlled by RKQS, which |
|---|
| 10 |
is called by RKCK45 (the mexfile) to perform the step [x1...x2] with |
|---|
| 11 |
4/5 order error control. RKCK45 controls the total integration from |
|---|
| 12 |
starting time to ending time. |
|---|
| 13 |
|
|---|
| 14 |
RKCK45 -->> |
|---|
| 15 |
RKQS -->> |
|---|
| 16 |
RKCK |
|---|
| 17 |
|
|---|
| 18 |
The following ANSI-C prototype describes how the rkqs wrapper RKCK45 |
|---|
| 19 |
is called, as suggested by NRC2.0. This is what this MEX is emulating, |
|---|
| 20 |
and will look like the MATLAB call to MATLAB's ode45, with some additional |
|---|
| 21 |
RHS arguments. We pass in the name of the m-file that |
|---|
| 22 |
contains the slopes (derivatives), just like in MATLAB's ode23(45) routines. |
|---|
| 23 |
|
|---|
| 24 |
void odeint(char *derivs_name name of m-file containing slopes |
|---|
| 25 |
double ystart[], initial conditions |
|---|
| 26 |
double x1, starting integration time |
|---|
| 27 |
double x2, ending integration time |
|---|
| 28 |
double eps, acceptable tolerance |
|---|
| 29 |
double h1, trial (hopeful) stepsize |
|---|
| 30 |
double hmin, minimum accetable stepsize |
|---|
| 31 |
int iret) number of wanted outputs |
|---|
| 32 |
|
|---|
| 33 |
TEST CALL: |
|---|
| 34 |
dname='dydx'; |
|---|
| 35 |
ic=[0 0 0 .00131]; |
|---|
| 36 |
t0=0.; |
|---|
| 37 |
tf=17.33*3600/2; |
|---|
| 38 |
err=1e-7; |
|---|
| 39 |
htry=900; |
|---|
| 40 |
hmin=1e-8; |
|---|
| 41 |
iret=1; |
|---|
| 42 |
[a,b,stat]=rkck45(dname,ic,t0,tf,err,htry,hmin,iret); |
|---|
| 43 |
|
|---|
| 44 |
dydx.m: |
|---|
| 45 |
function dydx = derivs(t,y) |
|---|
| 46 |
f=1e-4; |
|---|
| 47 |
dydx(1) = y(2); |
|---|
| 48 |
dydx(2) = f*y(4); |
|---|
| 49 |
dydx(3) = y(4); |
|---|
| 50 |
dydx(4) = -f*y(2); |
|---|
| 51 |
return |
|---|
| 52 |
|
|---|
| 53 |
a: 3.1194e+04 |
|---|
| 54 |
b: 26.1968 0.0000 0.2907 -0.0013 |
|---|
| 55 |
stat: |
|---|
| 56 |
1.0e+04 * {0.0032 |
|---|
| 57 |
0.0002 |
|---|
| 58 |
0.0001 |
|---|
| 59 |
3.1194 |
|---|
| 60 |
0 |
|---|
| 61 |
0 |
|---|
| 62 |
0 |
|---|
| 63 |
0 |
|---|
| 64 |
0 |
|---|
| 65 |
0.0004} |
|---|
| 66 |
|
|---|
| 67 |
|
|---|
| 68 |
|
|---|
| 69 |
Written by: |
|---|
| 70 |
BRIAN BLANTON |
|---|
| 71 |
UNIVERSITY OF NORTH CAROLINA AT CHAPEL HILL |
|---|
| 72 |
OCEAN PROCESSES NUMERICAL MODELING LABORATORY |
|---|
| 73 |
SPRING 97 |
|---|
| 74 |
|
|---|
| 75 |
*/ |
|---|
| 76 |
|
|---|
| 77 |
#include <math.h> |
|---|
| 78 |
#include <stdio.h> |
|---|
| 79 |
#include "mex.h" |
|---|
| 80 |
#include "opnml_mex5_allocs.c" |
|---|
| 81 |
#define SIGN(a,b) ((b) >= 0.0 ? fabs(a) : -fabs(a)) |
|---|
| 82 |
static double dmaxarg1,dmaxarg2; |
|---|
| 83 |
#define DMAX(a,b) (dmaxarg1=(a),dmaxarg2=(b),(dmaxarg1) > (dmaxarg2) ?\ |
|---|
| 84 |
(dmaxarg1) : (dmaxarg2)) |
|---|
| 85 |
#define MAXSTP 10000 |
|---|
| 86 |
#define TINY 1.0e-30 |
|---|
| 87 |
#define SAFETY 0.9 |
|---|
| 88 |
#define PGROW -0.2 |
|---|
| 89 |
#define PSHRNK -0.25 |
|---|
| 90 |
#define ERRCON 1.89e-4 |
|---|
| 91 |
|
|---|
| 92 |
/************************************************************ |
|---|
| 93 |
|
|---|
| 94 |
# ####### |
|---|
| 95 |
##### # # #### # # # # # |
|---|
| 96 |
# # # # # # # # # # # |
|---|
| 97 |
# # #### # #### # # ###### |
|---|
| 98 |
##### # # # # # ####### # |
|---|
| 99 |
# # # # # # # # # # # |
|---|
| 100 |
# # # # #### # # # ##### |
|---|
| 101 |
|
|---|
| 102 |
#### ## ##### ###### # # ## # # |
|---|
| 103 |
# # # # # # # # # # # # |
|---|
| 104 |
# # # # ##### # # # # # |
|---|
| 105 |
# ### ###### # # # ## # ###### # |
|---|
| 106 |
# # # # # # ## ## # # # |
|---|
| 107 |
#### # # # ###### # # # # # |
|---|
| 108 |
|
|---|
| 109 |
************************************************************/ |
|---|
| 110 |
|
|---|
| 111 |
/* PROTOTYPES */ |
|---|
| 112 |
void rkqs(double [], |
|---|
| 113 |
double [], |
|---|
| 114 |
int, |
|---|
| 115 |
double *, |
|---|
| 116 |
double, |
|---|
| 117 |
double, |
|---|
| 118 |
double [], |
|---|
| 119 |
double *, |
|---|
| 120 |
double *, |
|---|
| 121 |
char *); |
|---|
| 122 |
void rkck(double [], |
|---|
| 123 |
double [], |
|---|
| 124 |
int, |
|---|
| 125 |
double, |
|---|
| 126 |
double, |
|---|
| 127 |
double [], |
|---|
| 128 |
double [], |
|---|
| 129 |
char *); |
|---|
| 130 |
void help_fun(); |
|---|
| 131 |
void examp_fun(); |
|---|
| 132 |
void slopes_fun(); |
|---|
| 133 |
void stat_fun(); |
|---|
| 134 |
void disp_derivs(char *); |
|---|
| 135 |
|
|---|
| 136 |
/* Gateway/Entry */ |
|---|
| 137 |
void mexFunction(int nlhs, |
|---|
| 138 |
mxArray *plhs[], |
|---|
| 139 |
int nrhs, |
|---|
| 140 |
const mxArray *prhs[]) |
|---|
| 141 |
{ |
|---|
| 142 |
|
|---|
| 143 |
|
|---|
| 144 |
/* Declare RHS pointers; doubles for incoming from MATLAB, */ |
|---|
| 145 |
double *dx1,*dx2,*deps,*dh1,*dhmin,*dname; |
|---|
| 146 |
double *dkmax,*diret; |
|---|
| 147 |
int nok,nbad,nvar,kmax; |
|---|
| 148 |
|
|---|
| 149 |
/* Declare LHS pointers */ |
|---|
| 150 |
double *xpvec,*ypvec; |
|---|
| 151 |
double *stat; |
|---|
| 152 |
|
|---|
| 153 |
/* Local declarations */ |
|---|
| 154 |
double *xp,**yp; |
|---|
| 155 |
int i,idx,j,kount,m,n,nstp; |
|---|
| 156 |
double *ystart,x1,x2,eps,h1,hmin; |
|---|
| 157 |
double xsav,x,hnext,hdid,h,dxsav; |
|---|
| 158 |
double *yscal,*y,*dydx; |
|---|
| 159 |
char *slopes; |
|---|
| 160 |
mxArray *slp_rhs[2],*slp_lhs[1]; |
|---|
| 161 |
int errflag; |
|---|
| 162 |
|
|---|
| 163 |
/* ---- check I/O arguments ----------------------------------------- */ |
|---|
| 164 |
if (nrhs ==8){ |
|---|
| 165 |
if (nlhs !=2 & nlhs != 3) |
|---|
| 166 |
mexErrMsgTxt("rkck45 requires 2|3 output arguments."); |
|---|
| 167 |
|
|---|
| 168 |
/* First Matrix must contain a string */ |
|---|
| 169 |
if (!mxIsChar(prhs[0])) |
|---|
| 170 |
mexErrMsgTxt("RKCK45 must have string as first argument. HELP RKCK45"); |
|---|
| 171 |
|
|---|
| 172 |
/* Dereference first incoming (RHS) pointer */ |
|---|
| 173 |
/* Catch special calling options, like "help", "example", in first Matrix */ |
|---|
| 174 |
dname=mxGetPr(prhs[0]); |
|---|
| 175 |
m=mxGetN(prhs[0])+1; |
|---|
| 176 |
slopes=mxCalloc(m,sizeof(char)); |
|---|
| 177 |
mxGetString(prhs[0],slopes,m); |
|---|
| 178 |
/* Else, assume a valid mfile containing derivatives */ |
|---|
| 179 |
|
|---|
| 180 |
|
|---|
| 181 |
|
|---|
| 182 |
/* Dereference remaining incoming (RHS) pointers */ |
|---|
| 183 |
ystart=mxGetPr(prhs[1])-1; /* notice pointer decrement by 1 */ |
|---|
| 184 |
/* check dimensions of initial condition vector */ |
|---|
| 185 |
m=mxGetM(prhs[1]); |
|---|
| 186 |
n=mxGetN(prhs[1]); |
|---|
| 187 |
if(m*n != m & m*n != n ) |
|---|
| 188 |
mexErrMsgTxt("rkck45 requires initial conditions in vector form."); |
|---|
| 189 |
nvar=m*n; |
|---|
| 190 |
|
|---|
| 191 |
dx1=mxGetPr(prhs[2]);x1=dx1[0]; |
|---|
| 192 |
dx2=mxGetPr(prhs[3]);x2=dx2[0]; |
|---|
| 193 |
deps=mxGetPr(prhs[4]);eps=deps[0]; |
|---|
| 194 |
dh1=mxGetPr(prhs[5]);h1=dh1[0]; |
|---|
| 195 |
dhmin=mxGetPr(prhs[6]);hmin=dhmin[0]; |
|---|
| 196 |
diret=mxGetPr(prhs[7]);kmax=(int)diret[0]; |
|---|
| 197 |
if(kmax)dxsav=(x2-x1)/diret[0]; |
|---|
| 198 |
else dxsav=1; |
|---|
| 199 |
|
|---|
| 200 |
/* Local Allocations */ |
|---|
| 201 |
yscal =(double *) mxDvector(1,nvar); |
|---|
| 202 |
y =(double *) mxDvector(1,nvar); |
|---|
| 203 |
dydx =(double *) mxDvector(1,nvar); |
|---|
| 204 |
xp =(double *) mxDvector(1,MAXSTP); |
|---|
| 205 |
yp =(double **) mxDmatrix(1,MAXSTP,1,nvar); |
|---|
| 206 |
stat =(double *) mxDvector(0,4); |
|---|
| 207 |
|
|---|
| 208 |
x=x1; |
|---|
| 209 |
h=SIGN(h1,x2-x1); |
|---|
| 210 |
nok = nbad = kount = 0; |
|---|
| 211 |
for (i=1;i<=nvar;i++) y[i]=ystart[i]; |
|---|
| 212 |
if (kmax > 0) xsav=x-dxsav*2.0; /* Assure storage of first step */ |
|---|
| 213 |
mexSetTrapFlag(1); /* Disable exit-on-error behavior */ |
|---|
| 214 |
|
|---|
| 215 |
for (nstp=1;nstp<=MAXSTP;nstp++) { /* No more that MAXSTP steps */ |
|---|
| 216 |
|
|---|
| 217 |
/* Allocate Matrix structs for callback to MATLAB for derivatives */ |
|---|
| 218 |
slp_rhs[0]=mxCreateDoubleMatrix(1,1,mxREAL); |
|---|
| 219 |
slp_rhs[1]=mxCreateDoubleMatrix(1,nvar,mxREAL); |
|---|
| 220 |
slp_lhs[0]=mxCreateDoubleMatrix(1,nvar,mxREAL); |
|---|
| 221 |
mxSetPr(slp_rhs[0],&x); |
|---|
| 222 |
mxSetPr(slp_rhs[1],y+1); |
|---|
| 223 |
if (mexCallMATLAB(1,slp_lhs,2,slp_rhs,slopes)) /* Get starting slopes */ |
|---|
| 224 |
mexErrMsgTxt("call to mexCallMATLAB failed."); |
|---|
| 225 |
dydx=mxGetPr(slp_lhs[0])-1; |
|---|
| 226 |
|
|---|
| 227 |
/* Scaling used for accuracy monitoring */ |
|---|
| 228 |
for (i=1;i<=nvar;i++) yscal[i]=fabs(y[i])+fabs(dydx[i]*h)+TINY; |
|---|
| 229 |
if (kmax > 0 && kount < kmax-1 && fabs(x-xsav) > fabs(dxsav)) { |
|---|
| 230 |
xp[++kount]=x; |
|---|
| 231 |
for (i=1;i<=nvar;i++) yp[kount][i]=y[i]; |
|---|
| 232 |
xsav=x; |
|---|
| 233 |
} |
|---|
| 234 |
if ((x+h-x2)*(x+h-x1) > 0.0) h=x2-x; |
|---|
| 235 |
(*rkqs)(y,dydx,nvar,&x,h,eps,yscal,&hdid,&hnext,slopes); |
|---|
| 236 |
if (hdid == h) ++(nok); else ++(nbad); |
|---|
| 237 |
if ((x-x2)*(x2-x1) >= 0.0) { |
|---|
| 238 |
for (i=1;i<=nvar;i++) ystart[i]=y[i]; |
|---|
| 239 |
if (kmax) { |
|---|
| 240 |
xp[++kount]=x; |
|---|
| 241 |
for (i=1;i<=nvar;i++) yp[kount][i]=y[i]; |
|---|
| 242 |
|
|---|
| 243 |
/* This is where we need to store the intermediate |
|---|
| 244 |
results into the outgoing Matrices */ |
|---|
| 245 |
/* Create output matrices */ |
|---|
| 246 |
plhs[0]=mxCreateDoubleMatrix(kount,1,mxREAL); /* xp -> Matrix* plhs[0] */ |
|---|
| 247 |
plhs[1]=mxCreateDoubleMatrix(kount,nvar,mxREAL); /* yp -> Matrix* plhs[1] */ |
|---|
| 248 |
if(nlhs==3){ |
|---|
| 249 |
plhs[2]=mxCreateDoubleMatrix(5,1,mxREAL); /* status -> Matrix* plhs[1] */ |
|---|
| 250 |
stat[0]=nok; |
|---|
| 251 |
stat[1]=nbad; |
|---|
| 252 |
stat[2]=kount; |
|---|
| 253 |
stat[3]=dxsav; |
|---|
| 254 |
stat[4]=nvar; |
|---|
| 255 |
mxSetPr(plhs[2],stat); |
|---|
| 256 |
} |
|---|
| 257 |
xpvec=(double *) mxDvector(0,kount-1); |
|---|
| 258 |
ypvec=(double *) mxDvector(0,kount*nvar-1); |
|---|
| 259 |
|
|---|
| 260 |
for (i=1;i<=kount;i++)xpvec[i-1]=xp[i]; |
|---|
| 261 |
|
|---|
| 262 |
for (j=1;j<=nvar;j++){ |
|---|
| 263 |
for (i=1;i<=kount;i++){ |
|---|
| 264 |
idx=(j-1)*kount+i-1; |
|---|
| 265 |
ypvec[idx]=yp[i][j]; |
|---|
| 266 |
} |
|---|
| 267 |
} |
|---|
| 268 |
mxSetPr(plhs[0],xpvec); |
|---|
| 269 |
mxSetPr(plhs[1],ypvec); |
|---|
| 270 |
} |
|---|
| 271 |
else { |
|---|
| 272 |
plhs[0]=mxCreateDoubleMatrix(1,1,mxREAL); /* xp -> Matrix* plhs[0] */ |
|---|
| 273 |
plhs[1]=mxCreateDoubleMatrix(1,nvar,mxREAL); /* yp -> Matrix* plhs[1] */ |
|---|
| 274 |
xpvec=(double *)mxDvector(0,0); |
|---|
| 275 |
xpvec[0]=x; |
|---|
| 276 |
ypvec=(double *)mxDvector(0,nvar-1); |
|---|
| 277 |
for (j=1;j<=nvar;j++)ypvec[j-1]=y[j]; |
|---|
| 278 |
mxSetPr(plhs[0],xpvec); |
|---|
| 279 |
mxSetPr(plhs[1],ypvec); |
|---|
| 280 |
if(nlhs==3){ |
|---|
| 281 |
plhs[2]=mxCreateDoubleMatrix(5,1,mxREAL); /* status -> Matrix* plhs[1] */ |
|---|
| 282 |
stat[0]=nok; |
|---|
| 283 |
stat[1]=nbad; |
|---|
| 284 |
stat[4]=nvar; |
|---|
| 285 |
mxSetPr(plhs[2],stat); |
|---|
| 286 |
} |
|---|
| 287 |
} |
|---|
| 288 |
return; |
|---|
| 289 |
} |
|---|
| 290 |
if (fabs(hnext) <= hmin) { |
|---|
| 291 |
fprintf(stderr,"Next Stepsize (%e) exceeds Minimum (%e)\n",hnext,hmin); |
|---|
| 292 |
mexErrMsgTxt("Step size too small in RKCK45"); |
|---|
| 293 |
} |
|---|
| 294 |
h=hnext; |
|---|
| 295 |
} |
|---|
| 296 |
mexErrMsgTxt("Too many steps in routine odeint"); |
|---|
| 297 |
} |
|---|
| 298 |
else if(nrhs==1) |
|---|
| 299 |
{ |
|---|
| 300 |
/* First Matrix must contain a string */ |
|---|
| 301 |
if (!mxIsChar(prhs[0])) |
|---|
| 302 |
mexErrMsgTxt("RKCK45 must have string as first argument. HELP RKCK45"); |
|---|
| 303 |
/* Get option from prhs[0] and process */ |
|---|
| 304 |
dname=mxGetPr(prhs[0]); |
|---|
| 305 |
m=mxGetN(prhs[0])+1; |
|---|
| 306 |
slopes=mxCalloc(m,sizeof(char)); |
|---|
| 307 |
mxGetString(prhs[0],slopes,m); |
|---|
| 308 |
if(strcmp(slopes,"help")==0)help_fun(); |
|---|
| 309 |
else if(strcmp(slopes,"example")==0)examp_fun(); |
|---|
| 310 |
else if(strcmp(slopes,"slopes")==0)slopes_fun(); |
|---|
| 311 |
else if(strcmp(slopes,"stat")==0)stat_fun(); |
|---|
| 312 |
else |
|---|
| 313 |
{ |
|---|
| 314 |
/* Try to determine if the value of "slopes" points to a valid mfile */ |
|---|
| 315 |
mxArray *xrhs[1],*xlhs[1]; |
|---|
| 316 |
double *derr; |
|---|
| 317 |
int err; |
|---|
| 318 |
puts(slopes); |
|---|
| 319 |
xrhs[0]=mxCreateString(slopes); |
|---|
| 320 |
xlhs[0]=mxCreateDoubleMatrix(1,1,mxREAL); |
|---|
| 321 |
mexCallMATLAB(1,xlhs,1,xrhs,"exist"); /* Probe Matlab */ |
|---|
| 322 |
derr=mxGetPr(xlhs[0]); |
|---|
| 323 |
err=(int)derr[0]; |
|---|
| 324 |
if(err==2) |
|---|
| 325 |
disp_derivs(slopes); |
|---|
| 326 |
else |
|---|
| 327 |
mexErrMsgTxt("Unrecognized first argument to RKCK45"); |
|---|
| 328 |
} |
|---|
| 329 |
|
|---|
| 330 |
} |
|---|
| 331 |
else |
|---|
| 332 |
mexErrMsgTxt("RKCK45 requires 1|8 input (RHS) arguments"); |
|---|
| 333 |
} |
|---|
| 334 |
|
|---|
| 335 |
void rkqs(double y[], double dydx[], int n, |
|---|
| 336 |
double *x, double htry, double eps, |
|---|
| 337 |
double yscal[], double *hdid, |
|---|
| 338 |
double *hnext,char *slopes) |
|---|
| 339 |
{ |
|---|
| 340 |
int i; |
|---|
| 341 |
double errmax,h,xnew,*yerr,*ytemp; |
|---|
| 342 |
yerr=mxDvector(1,n); |
|---|
| 343 |
ytemp=mxDvector(1,n); |
|---|
| 344 |
h=htry; |
|---|
| 345 |
for (;;) { |
|---|
| 346 |
rkck(y,dydx,n,*x,h,ytemp,yerr,slopes); |
|---|
| 347 |
errmax=0.0; |
|---|
| 348 |
for (i=1;i<=n;i++) errmax=DMAX(errmax,fabs(yerr[i]/yscal[i])); |
|---|
| 349 |
errmax /= eps; |
|---|
| 350 |
if (errmax > 1.0) { |
|---|
| 351 |
h=SAFETY*h*pow(errmax,PSHRNK); |
|---|
| 352 |
if (h < 0.1*h) h *= 0.1; |
|---|
| 353 |
xnew=(*x)+h; |
|---|
| 354 |
if (xnew == *x) mexErrMsgTxt("stepsize underflow in rkqs"); |
|---|
| 355 |
continue; |
|---|
| 356 |
} else { |
|---|
| 357 |
if (errmax > ERRCON) *hnext=SAFETY*h*pow(errmax,PGROW); |
|---|
| 358 |
else *hnext=5.0*h; |
|---|
| 359 |
*x += (*hdid=h); |
|---|
| 360 |
for (i=1;i<=n;i++) y[i]=ytemp[i]; |
|---|
| 361 |
break; |
|---|
| 362 |
} |
|---|
| 363 |
} |
|---|
| 364 |
} |
|---|
| 365 |
|
|---|
| 366 |
void rkck(double y[], double dydx[], int n, |
|---|
| 367 |
double x, double h, double yout[], |
|---|
| 368 |
double yerr[],char *slopes) |
|---|
| 369 |
{ |
|---|
| 370 |
int i; |
|---|
| 371 |
static double a2=0.2,a3=0.3,a4=0.6,a5=1.0,a6=0.875,b21=0.2, |
|---|
| 372 |
b31=3.0/40.0,b32=9.0/40.0,b41=0.3,b42 = -0.9,b43=1.2, |
|---|
| 373 |
b51 = -11.0/54.0, b52=2.5,b53 = -70.0/27.0,b54=35.0/27.0, |
|---|
| 374 |
b61=1631.0/55296.0,b62=175.0/512.0,b63=575.0/13824.0, |
|---|
| 375 |
b64=44275.0/110592.0,b65=253.0/4096.0,c1=37.0/378.0, |
|---|
| 376 |
c3=250.0/621.0,c4=125.0/594.0,c6=512.0/1771.0, |
|---|
| 377 |
dc5 = -277.0/14336.0; |
|---|
| 378 |
double dc1=c1-2825.0/27648.0,dc3=c3-18575.0/48384.0, |
|---|
| 379 |
dc4=c4-13525.0/55296.0,dc6=c6-0.25; |
|---|
| 380 |
double *ak2,*ak3,*ak4,*ak5,*ak6,*ytemp; |
|---|
| 381 |
mxArray *slp_rhs[2],*slp_lhs[1]; |
|---|
| 382 |
double x2,x3,x4,x5,x6; |
|---|
| 383 |
|
|---|
| 384 |
x2=x+a2*h; |
|---|
| 385 |
x3=x+a3*h; |
|---|
| 386 |
x4=x+a4*h; |
|---|
| 387 |
x5=x+a5*h; |
|---|
| 388 |
x6=x+a6*h; |
|---|
| 389 |
|
|---|
| 390 |
slp_rhs[0]=mxCreateDoubleMatrix(1,1,mxREAL); |
|---|
| 391 |
slp_rhs[1]=mxCreateDoubleMatrix(1,n,mxREAL); |
|---|
| 392 |
slp_lhs[0]=mxCreateDoubleMatrix(1,n,mxREAL); |
|---|
| 393 |
ak2=mxDvector(1,n); |
|---|
| 394 |
ak3=mxDvector(1,n); |
|---|
| 395 |
ak4=mxDvector(1,n); |
|---|
| 396 |
ak5=mxDvector(1,n); |
|---|
| 397 |
ak6=mxDvector(1,n); |
|---|
| 398 |
ytemp=mxDvector(1,n); |
|---|
| 399 |
mxSetPr(slp_rhs[1],ytemp+1); |
|---|
| 400 |
|
|---|
| 401 |
for (i=1;i<=n;i++) |
|---|
| 402 |
ytemp[i]=y[i]+b21*h*dydx[i]; |
|---|
| 403 |
|
|---|
| 404 |
mxSetPr(slp_rhs[0],&x2); |
|---|
| 405 |
if (mexCallMATLAB(1,slp_lhs,2,slp_rhs,slopes)) |
|---|
| 406 |
mexErrMsgTxt("call to mexCallMATLAB failed."); |
|---|
| 407 |
ak2=mxGetPr(slp_lhs[0])-1; |
|---|
| 408 |
for (i=1;i<=n;i++) |
|---|
| 409 |
ytemp[i]=y[i]+h*(b31*dydx[i]+b32*ak2[i]); |
|---|
| 410 |
|
|---|
| 411 |
mxSetPr(slp_rhs[0],&x3); |
|---|
| 412 |
if (mexCallMATLAB(1,slp_lhs,2,slp_rhs,slopes)) |
|---|
| 413 |
mexErrMsgTxt("call to mexCallMATLAB failed."); |
|---|
| 414 |
ak3=mxGetPr(slp_lhs[0])-1; |
|---|
| 415 |
for (i=1;i<=n;i++) |
|---|
| 416 |
ytemp[i]=y[i]+h*(b41*dydx[i]+b42*ak2[i]+b43*ak3[i]); |
|---|
| 417 |
|
|---|
| 418 |
mxSetPr(slp_rhs[0],&x4); |
|---|
| 419 |
if (mexCallMATLAB(1,slp_lhs,2,slp_rhs,slopes)) |
|---|
| 420 |
mexErrMsgTxt("call to mexCallMATLAB failed."); |
|---|
| 421 |
ak4=mxGetPr(slp_lhs[0])-1; |
|---|
| 422 |
for (i=1;i<=n;i++) |
|---|
| 423 |
ytemp[i]=y[i]+h*(b51*dydx[i]+b52*ak2[i]+b53*ak3[i]+b54*ak4[i]); |
|---|
| 424 |
|
|---|
| 425 |
mxSetPr(slp_rhs[0],&x5); |
|---|
| 426 |
if (mexCallMATLAB(1,slp_lhs,2,slp_rhs,slopes)) |
|---|
| 427 |
mexErrMsgTxt("call to mexCallMATLAB failed."); |
|---|
| 428 |
ak5=mxGetPr(slp_lhs[0])-1; |
|---|
| 429 |
for (i=1;i<=n;i++) |
|---|
| 430 |
ytemp[i]=y[i]+h*(b61*dydx[i]+b62*ak2[i]+b63*ak3[i]+b64*ak4[i]+b65*ak5[i]); |
|---|
| 431 |
|
|---|
| 432 |
mxSetPr(slp_rhs[0],&x6); |
|---|
| 433 |
if (mexCallMATLAB(1,slp_lhs,2,slp_rhs,slopes)) |
|---|
| 434 |
mexErrMsgTxt("call to mexCallMATLAB failed."); |
|---|
| 435 |
ak6=mxGetPr(slp_lhs[0])-1; |
|---|
| 436 |
for (i=1;i<=n;i++) |
|---|
| 437 |
yout[i]=y[i]+h*(c1*dydx[i]+c3*ak3[i]+c4*ak4[i]+c6*ak6[i]); |
|---|
| 438 |
|
|---|
| 439 |
for (i=1;i<=n;i++) |
|---|
| 440 |
yerr[i]=h*(dc1*dydx[i]+dc3*ak3[i]+dc4*ak4[i]+dc5*ak5[i]+dc6*ak6[i]); |
|---|
| 441 |
} |
|---|
| 442 |
|
|---|
| 443 |
void help_fun() |
|---|
| 444 |
{ |
|---|
| 445 |
fprintf(stderr,"\n"); |
|---|
| 446 |
fprintf(stderr," This is the HELP page for RKCK45\n"); |
|---|
| 447 |
fprintf(stderr," See Numerical Recipes in C, 2nd Ed.\n"); |
|---|
| 448 |
fprintf(stderr," Chapter 16, Sections 1 and 2.\n"); |
|---|
| 449 |
fprintf(stderr,"\n"); |
|---|
| 450 |
|
|---|
| 451 |
} |
|---|
| 452 |
|
|---|
| 453 |
void stat_fun() |
|---|
| 454 |
{ |
|---|
| 455 |
fprintf(stderr," This is the STAT page for RKCK45\n"); |
|---|
| 456 |
fprintf(stderr,"\n"); |
|---|
| 457 |
fprintf(stderr," The 'stat' return variable is optional.\n"); |
|---|
| 458 |
fprintf(stderr," If returned, it contains the following entries:\n"); |
|---|
| 459 |
fprintf(stderr," stat(1) -> Number of successful timesteps (nok)\n"); |
|---|
| 460 |
fprintf(stderr," stat(2) -> Number of unsuccessful timesteps (nbad)\n"); |
|---|
| 461 |
fprintf(stderr," stat(3) -> Number of intermediate results returned\n"); |
|---|
| 462 |
fprintf(stderr," stat(4) -> Output interval, not really useful.\n"); |
|---|
| 463 |
fprintf(stderr," stat(5) -> Number of equations\n"); |
|---|
| 464 |
fprintf(stderr,"\n"); |
|---|
| 465 |
} |
|---|
| 466 |
|
|---|
| 467 |
void examp_fun() |
|---|
| 468 |
{ |
|---|
| 469 |
fprintf(stderr," This is the EXAMPLE page for RKCK45\n"); |
|---|
| 470 |
fprintf(stderr,"\n"); |
|---|
| 471 |
fprintf(stderr," The following example is actually the integration of\n"); |
|---|
| 472 |
fprintf(stderr," the equations of motion describing the translation of\n"); |
|---|
| 473 |
fprintf(stderr," a particle on a rotating, frictionless plane. How about\n"); |
|---|
| 474 |
fprintf(stderr," a hockey puck struck with sufficient velocity (m/s) to\n"); |
|---|
| 475 |
fprintf(stderr," maintain an inertial circle within a hockey rink 26 m\n"); |
|---|
| 476 |
fprintf(stderr," wide in Savannah, Georgia.\n"); |
|---|
| 477 |
fprintf(stderr,"\n"); |
|---|
| 478 |
fprintf(stderr," Given the following derivatives mfile dydx.m:\n"); |
|---|
| 479 |
fprintf(stderr,"\n"); |
|---|
| 480 |
fprintf(stderr," function dydx = derivs(t,y)\n"); |
|---|
| 481 |
fprintf(stderr," f=7.7368e-05; /* Coriolis */\n"); |
|---|
| 482 |
fprintf(stderr," dydx(1) = y(2);\n"); |
|---|
| 483 |
fprintf(stderr," dydx(2) = f*y(4);\n"); |
|---|
| 484 |
fprintf(stderr," dydx(3) = y(4);\n"); |
|---|
| 485 |
fprintf(stderr," dydx(4) = -f*y(2);\n"); |
|---|
| 486 |
fprintf(stderr," return\n"); |
|---|
| 487 |
fprintf(stderr,"\n"); |
|---|
| 488 |
|
|---|
| 489 |
fprintf(stderr," and the following parameters:\n"); |
|---|
| 490 |
fprintf(stderr," dname='dydx'; /* Name of derivatives mfile */\n"); |
|---|
| 491 |
fprintf(stderr," ic=[0 0 0 .00101]; /* Xo,Uo,Yo,Vo */\n"); |
|---|
| 492 |
fprintf(stderr," t0=0.;\n"); |
|---|
| 493 |
fprintf(stderr," tf=22.5587*3600/2; /* This is .5*T_inertial */\n"); |
|---|
| 494 |
fprintf(stderr," err=1e-7;\n"); |
|---|
| 495 |
fprintf(stderr," htry=90; /* Desired stepsize */\n"); |
|---|
| 496 |
fprintf(stderr," hmin=1e-8; /* Minimum acceptable stepsize */\n"); |
|---|
| 497 |
fprintf(stderr," iret=25; /* Keep at most 25 intermediate results */\n"); |
|---|
| 498 |
fprintf(stderr,"\n"); |
|---|
| 499 |
|
|---|
| 500 |
fprintf(stderr," the following call:\n"); |
|---|
| 501 |
fprintf(stderr," [a,b,stat]=rkck45(dname,ic,t0,tf,err,htry,hmin,iret);\n"); |
|---|
| 502 |
fprintf(stderr,"\n"); |
|---|
| 503 |
|
|---|
| 504 |
fprintf(stderr," produces these results:\n"); |
|---|
| 505 |
|
|---|
| 506 |
fprintf(stderr," a= 4.0606e+04\n"); |
|---|
| 507 |
fprintf(stderr," b= 26.1090 0.0000 0.0002 -0.0010\n"); |
|---|
| 508 |
fprintf(stderr," stat= \n"); |
|---|
| 509 |
fprintf(stderr," 1.0e+04 * 0.0027\n"); |
|---|
| 510 |
fprintf(stderr," 0.0002\n"); |
|---|
| 511 |
fprintf(stderr," 0.0001\n"); |
|---|
| 512 |
fprintf(stderr," 4.0606\n"); |
|---|
| 513 |
fprintf(stderr," 0.0004\n"); |
|---|
| 514 |
fprintf(stderr,"\n"); |
|---|
| 515 |
} |
|---|
| 516 |
void slopes_fun() |
|---|
| 517 |
{ |
|---|
| 518 |
fprintf(stderr,"\n"); |
|---|
| 519 |
fprintf(stderr," This is the SLOPES help for RKCK45\n"); |
|---|
| 520 |
fprintf(stderr," Example derivatives mfile called dydx.m:\n"); |
|---|
| 521 |
fprintf(stderr,"\n"); |
|---|
| 522 |
fprintf(stderr," function dydx = derivs(t,y)\n"); |
|---|
| 523 |
fprintf(stderr," f=1e-4; \n"); |
|---|
| 524 |
fprintf(stderr," dydx(1) = y(2);\n"); |
|---|
| 525 |
fprintf(stderr," dydx(2) = f*y(4);\n"); |
|---|
| 526 |
fprintf(stderr," dydx(3) = y(4);\n"); |
|---|
| 527 |
fprintf(stderr," dydx(4) = -f*y(2);\n"); |
|---|
| 528 |
fprintf(stderr," return\n"); |
|---|
| 529 |
fprintf(stderr,"\n"); |
|---|
| 530 |
|
|---|
| 531 |
} |
|---|
| 532 |
void disp_derivs(char *slopes) |
|---|
| 533 |
{ |
|---|
| 534 |
mxArray *xrhs[1],*xlhs[1]; |
|---|
| 535 |
double *derr; |
|---|
| 536 |
int err; |
|---|
| 537 |
xrhs[0]=mxCreateString(slopes); |
|---|
| 538 |
xlhs[0]=mxCreateDoubleMatrix(1,1,mxREAL); |
|---|
| 539 |
fprintf(stderr,"\n"); |
|---|
| 540 |
fprintf(stderr," Your designated derivatives file looks like:\n"); |
|---|
| 541 |
|
|---|
| 542 |
mexCallMATLAB(0,xlhs,1,xrhs,"type"); /* Probe Matlab */ |
|---|
| 543 |
|
|---|
| 544 |
} |
|---|
| 545 |
#undef SAFETY |
|---|
| 546 |
#undef PGROW |
|---|
| 547 |
#undef PSHRNK |
|---|
| 548 |
#undef ERRCON |
|---|
