1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37
|
void runge_kutta_evl(unsigned int I3, double Dt, vector<double>* T, vector<double>* Y1, vector<double>* Y2, vector<double>* Y3, vector<double>* Y5, vector<double>* r_range, vector<double>* m_range,vector<double>* rho_range){
unsigned int pp_size=r_range.size();
double m_int_0=Linear_once(pp_size,r_range, m_range,(*Y2)[I3]);
double rho_int_0=Linear_once(pp_size,r_range, rho_range,(*Y2)[I3]);
double a0_y1 = Dt*dy1_dt((*Y1)[I3], (*Y2)[I3], (*Y3)[I3], (*Y5)[I3], m_int_0, rho_int_0);
double b0_y2=Dt*dy2_dt((*Y1)[I3]);
double c0_y3=Dt*dy3_dt((*Y1)[I3], (*Y2)[I3], (*Y5)[I3], rho_int_0);
double e0_y5=Dt*dy5_dt((*Y1)[I3], (*Y2)[I3], (*Y3)[I3], (*Y5)[I3], rho_int_0);
double m_int_1=Linear_once(pp_size,r_range, m_range,(*Y2)[I3]+0.5*b0_y2);
double rho_int_1=Linear_once(pp_size,r_range, rho_range,(*Y2)[I3]+0.5*b0_y2);
double a1_y1 = Dt*dy1_dt((*Y1)[I3]+0.5*a0_y1, (*Y2)[I3]+0.5*b0_y2, (*Y3)[I3]+0.5*c0_y3, (*Y5)[I3]+0.5*e0_y5, m_int_1, rho_int_1);
double b1_y2 = Dt*dy2_dt((*Y1)[I3]+0.5*a0_y1);
double c1_y3 = Dt*dy3_dt((*Y1)[I3]+0.5*a0_y1, (*Y2)[I3]+0.5*b0_y2, (*Y5)[I3]+0.5*e0_y5, rho_int_1);
double e1_y5 = Dt*dy5_dt((*Y1)[I3]+0.5*a0_y1, (*Y2)[I3]+0.5*b0_y2, (*Y3)[I3]+0.5*c0_y3, (*Y5)[I3]+0.5*e0_y5, rho_int_1);
double m_int_2=Linear_once(pp_size,r_range, m_range,(*Y2)[I3]+0.5*b1_y2);
double rho_int_2=Linear_once(pp_size,r_range, rho_range, (*Y2)[I3]+0.5*b1_y2);
double a2_y1 = Dt*dy1_dt((*Y1)[I3]+0.5*a1_y1, (*Y2)[I3]+0.5*b1_y2, (*Y3)[I3]+0.5*c1_y3, (*Y5)[I3]+0.5*e1_y5, m_int_2, rho_int_2);
double b2_y2 = Dt*dy2_dt((*Y1)[I3]+0.5*a1_y1);
double c2_y3 = Dt*dy3_dt((*Y1)[I3]+0.5*a1_y1, (*Y2)[I3]+0.5*b1_y2, (*Y5)[I3]+0.5*e1_y5, rho_int_2);
double e2_y5 = Dt*dy5_dt((*Y1)[I3]+0.5*a1_y1, (*Y2)[I3]+0.5*b1_y2, (*Y3)[I3]+0.5*c1_y3, (*Y5)[I3]+0.5*e1_y5, rho_int_2);
double m_int_3=Linear_once(pp_size, r_range, m_range, (*Y2)[I3]+b2_y2);
double rho_int_3=Linear_once(pp_size,r_range, rho_range, (*Y2)[I3]+b2_y2);
double a3_y1 = Dt*dy1_dt((*Y1)[I3]+a2_y1, (*Y2)[I3]+b2_y2, (*Y3)[I3]+c2_y3, (*Y5)[I3]+e2_y5, m_int_3, rho_int_3);
double b3_y2 = Dt*dy2_dt((*Y1)[I3]+a2_y1);
double c3_y3 = Dt*dy3_dt((*Y1)[I3]+a2_y1, (*Y2)[I3]+b2_y2, (*Y5)[I3]+e2_y5, rho_int_3);
double e3_y5 = Dt*dy5_dt((*Y1)[I3]+a2_y1, (*Y2)[I3]+b2_y2, (*Y3)[I3]+c2_y3, (*Y5)[I3]+e2_y5, rho_int_3);
*T.push_back((*T)[0] + Dt);
(*Y1).push_back((*Y1)[I3] + (a0_y1 + 2*a1_y1 + 2*a2_y1 + a3_y1)/6);
(*Y2).push_back((*Y2)[I3] + (b0_y2 + 2*b1_y2 + 2*b2_y2 + b3_y2)/6);
(*Y3).push_back((*Y3)[I3] + (c0_y3 + 2*c1_y3 + 2*c2_y3 + c3_y3)/6);
(*Y5).push_back((*Y5)[I3] + (e0_y5 + 2*e1_y5 + 2*e2_y5 + e3_y5)/6);
}
|