32 lines
575 B
C++
32 lines
575 B
C++
#include "nr.h"
|
|
|
|
void NR::mmid(Vec_I_DP &y, Vec_I_DP &dydx, const DP xs, const DP htot,
|
|
const int nstep, Vec_O_DP &yout,
|
|
void derivs(const DP, Vec_I_DP &, Vec_O_DP &))
|
|
{
|
|
int i,n;
|
|
DP x,swap,h2,h;
|
|
|
|
int nvar=y.size();
|
|
Vec_DP ym(nvar),yn(nvar);
|
|
h=htot/nstep;
|
|
for (i=0;i<nvar;i++) {
|
|
ym[i]=y[i];
|
|
yn[i]=y[i]+h*dydx[i];
|
|
}
|
|
x=xs+h;
|
|
derivs(x,yn,yout);
|
|
h2=2.0*h;
|
|
for (n=1;n<nstep;n++) {
|
|
for (i=0;i<nvar;i++) {
|
|
swap=ym[i]+h2*yout[i];
|
|
ym[i]=yn[i];
|
|
yn[i]=swap;
|
|
}
|
|
x += h;
|
|
derivs(x,yn,yout);
|
|
}
|
|
for (i=0;i<nvar;i++)
|
|
yout[i]=0.5*(ym[i]+yn[i]+h*yout[i]);
|
|
}
|