Can anyone help please? I would love a quick reply, I need this asap.

Thanks.

I know how to scan an array and test a variable eg: if min > array[i] set min to array[i], but i do not know how to get energy(y) into an array like array[i].

here's the program:

#include <iostream> #include <cmath> #include <algorithm> using namespace std; void rhs (double rho[], double f[]); void step_RK2(double h, double yin[], double yout[]); double energy (double y[]); main() { double rhox,vx,tau; double rhoy,vy; double h,hnew; double E0; double y[4],yh[4],y2h[4]; int i; double safety=0.9,max_increase=2,max_decrease=5; double ratio,eps,eps0,atol,rtol; int accept,p,two_to_p; h=0.01; atol=rtol=1e-3; tau=0; rhox=0;rhoy=1; vy=0.4; vx=0.4; /* the following lines initialise array y; they also serve as a translation table between vector elements and physical variables */ y[0]=rhox; y[1]=rhoy; y[2]=vx; y[3]=vy; E0=energy(y); cout << tau << "\t" << rhox << " \t" << rhoy << " \t" << vx << " \t" << vy << " \t" << E0 << " \t" << 0 << " \t" << h << "\n" ; while (tau < 300) { /* two smaller steps */ step_RK2(h,y ,yh); step_RK2(h,yh,yh); /* one double step */ step_RK2(2*h,y ,y2h); /* error estimate for 2nd-order method*/ p=2;two_to_p=4; ratio=0; for(i=0;i<4;i++) { eps=abs(y2h[i] - yh[i])/(two_to_p - 1); eps0=atol+abs(yh[i])*rtol; ratio=max(ratio,eps/eps0); } /* step size correction */ hnew=h/exp( log(ratio)/(p+1) ); /* accept/reject the step */ if(ratio<1) { if(accept==1) h=min(max_increase*h,safety*hnew); for(i=0;i<4;i++) y[i]=(two_to_p*yh[i] - y2h[i])/(two_to_p - 1); accept=1; tau += 2*h; } else { accept=0; h=max(h/max_decrease,safety*hnew); }; cout << tau << "\t" << y[0] << " \t" << y[1] << " \t" << y[2] << " \t" << y[3] << " \t" << energy(y) << " \t" << energy(y)-E0 << " \t" << h << "\n" ; } ; } ; void rhs (double y[], double force[]) { double rho_abs,rho3; int i; rho_abs=sqrt(y[0]*y[0]+y[1]*y[1]); rho3=rho_abs*rho_abs*rho_abs; force[0] = y[2]; /* vx */ force[1] = y[3]; /* vx */ force[2] = -y[0]/rho3; /* -rhox/rho3 */ force[3] = -y[1]/rho3; /* -rhoy/rho3 */ }; void step_RK2(double h, double yin[], double yout[]) { double ytemp[4],k1[4],k2[4]; int i; /* evaluate k1 */ rhs(yin,k1); /* intermediate step */ for(i=0;i<4;i++) ytemp[i] = yin[i] + h*k1[i]; /* evaluate k2 */ rhs(ytemp,k2); /* update y */ for(i=0;i<4;i++) yout[i] = yin[i] + h*(k1[i]+k2[i])/2; }; double energy (double y[]) { return( (y[2]*y[2]+y[3]*y[3])/2 - 1/sqrt(y[0]*y[0]+y[1]*y[1]) ); }

This post has been edited by **Jayman**: 26 October 2010 - 06:47 PM