Problema de aplicación: métodos Runge-Kutta 1, 3 y 5.

Este programa resuelve numéricamente, mediante los métodos RK 1, RK2 y RK3 el problema de calcular la rapidez en un instante de tiempo “t” de una partícula relativista que pierde masa a una tasa constante, hasta llegar a una masa mínima. El medio ejerce dos tipos de fricción sobre la partícula: una proporcional a su masa y otra proporcional a su rapidez. El modelo físico a resolver es:

d(m[t] g[t] v[t])/dt = – ( k m[t] + b v[t] ),

donde

m[t] = mo – at, y g[t] es el factor gama relativista.

Como resultados, se obtiene que si la tasa de pérdida de masa multiplicada por el factor gama es igual al coeficiente de fricción con el medio (dependiente de la rapidez), entonces la partícula no gana rapidez como resultado de perder masa, experimentando desaceleración y frenándose con el tiempo. Si no hay fricción de ningún tipo, la partícula se acelera hasta alcanzar la velocidad de la luz. Este fenómeno es consecuencia de la ley de conservación del momento lineal.

Como era de esperarse, la rapidez de la partícula disminuye una vez que ha alcanzado su masa mínima.

Los métodos RK3 y RK5 arrojaron resultados cercanos entre sí, mientras que el método RK1 arrojó resultados cercanos a los de los otros dos métodos sólo en los extremos, aunque la forma de la curva se aproxima. No se aprecian discontinuidades en las curvas.

Como observación, debe señalarse que es necesario ajustar el tamaño de paso a un valor pequeño, ya que se corre el riesgo de obtener una masa negativa en pocas iteraciones, lo cual arruina los resultados. El tamaño de paso debe ser tal que, con una tasa de pérdida de masa dada, m[t] alcance exactamente el valor de la masa mínima definida en el código.

Una gráfica con las curvas obtenidas puede verse arriba y hasta abajo de este post.

================================== Código ==================================

#include <iostream>
#include <cmath>
#include <fstream>

using namespace std;

int main(){
// condiciones iniciales
double v1 = 2.9*pow(10,8);
double v3 = v1;
double v5 = v3;
double t = 0;
int i = 0;

// definicion de coeficientes
double k1, k2, k3, k4, k5, k6, kt, kr1, kr2, kr3, krt;

// definicion de parametros
double c = 3*pow(10,8);  // velocidad de la luz
double a = 1;            // tasa de perdida de masa kg/s
double m = 50;            // masa inicial kg
double mmin = 0.1;       // masa minima kg
double k = 1;            // coeficiente de friccion masica
double b = 0.5;          // coeficiente de friccion por velocidad
double dt = 0.1;         // tamaño de paso

// para escribir output
ofstream file;
file.open(“rk-caja.txt”);

cout << “Metodos RK:”<<endl;
cout << “it \t t \t m \t v(rk1) \t v(rk3) \t v(rk5)” << endl << endl;

file << “Metodos RK:”<<endl;
file << “it \t t \t m \t v(rk1) \t v(rk3) \t v(rk5)” << endl << endl;

//————— región para cuando la masa es mayor a la masa mínima ——————–//
while(m – a*t > mmin){

// ————————————— RK 1 ——————————————————————-//
k = pow(1 – pow(v1/c,2),3/2)*((v1/(m-a*t))*(a*pow(1-pow(v1/c,2),-1/2) – b) – k)*dt;

// ————————————— RK 3 ——————————————————————-//
kr1 = pow(1 – pow(v3/c,2),3/2)*((v3/(m-a*t))*(a*pow(1-pow(v3/c,2),-1/2) – b) – k)*dt;
kr2 = pow(1 – pow((v3+k1/3)/c,2),3/2)*(((v3+k1/3)/(m-a*(t+dt/3)))*(a*pow(1-pow((v3+k1/3)/c,2),-1/2) – b) – k)*dt;
kr3 = pow(1 – pow((v3+2*k2/3)/c,2),3/2)*(((v3+2*k2/3)/(m-a*(t+2*dt/3)))*(a*pow(1-pow((v3+2*k2/3)/c,2),-1/2) – b) – k)*dt;
krt = (kr1 + 3*kr3)/4;

// ————————————— RK 5 ——————————————————————-//
k1 = pow(1 – pow(v5/c,2),3/2)*((v5/(m-a*t))*(a*pow(1-pow(v5/c,2),-1/2) – b) – k)*dt;
k2 = pow(1 – pow((v5+k1)/c,2),3/2)*(((v5+k1)/(m-a*(t+dt)))*(a*pow(1-pow((v5+k1)/c,2),-1/2) – b) – k)*dt;
k3 = pow(1 – pow((v5+(k1+k2)/2)/c,2),3/2)*(((v5+(k1+k2)/2)/(m-a*(t+dt)))*(a*pow(1-pow((v5+(k1+k2)/2)/c,2),-1/2) – b) – k)*dt;
k4 = pow(1 – pow((v5+(14*k1+5*k2-3*k3)/64)/c,2),3/2)*(((v5+(14*k1+5*k2-3*k3)/64)/(m-a*(t+dt/4)))*(a*pow(1-pow((v5+(14*k1+5*k2-3*k3)/64)/c,2),-1/2) – b) – k)*dt;
k5 = pow(1 – pow((v5+(-12*k1-12*k2+8*k3+64*k4)/96)/c,2),3/2)*(((v5+(-12*k1-12*k2+8*k3+64*k4)/96)/(m-a*(t+dt/2)))*(a*pow(1-pow((v5+(-12*k1-12*k2+8*k3+64*k4)/96)/c,2),-1/2) – b) – k)*dt;
k6 = pow(1 – pow((v5+(-9*k2 +5*k3 +16*k4 +36*k5)/64)/c,2),3/2)*(((v5+(-9*k2 +5*k3 +16*k4 +36*k5)/64)/(m-a*(t+3*dt/4)))*(a*pow(1-pow((v5+(-9*k2 +5*k3 +16*k4 +36*k5)/64)/c,2),-1/2) – b) – k)*dt;
kt = (7*k1 + 7*k3 +32*k4 +12*k5 +32*k6)/90;

cout << i << “\t” << t << “\t” << m – a*t << “\t” << v1 << “\t” << v3 << “\t” << v5 << “\t” <<endl;
file << i << “\t” << t << “\t” << m – a*t << “\t” << v1 << “\t” << v3 << “\t” << v5 << “\t” << endl;

v1 += k;
v3 += krt;
v5 += kt;

t += dt;
i ++;
}
//—————- región para cuando la partícula ya no puede perder masa ——————-//
while(v5 > 0){

// ————————————— RK 1 ——————————————————————-//
k = – pow(1 – pow(v1/c,2),3/2)*((v1/mmin)*b + k)*dt;

// ————————————— RK 3 ——————————————————————-//
kr1 = – pow(1 – pow(v3/c,2),3/2)*((v3/mmin)*b + k)*dt;
kr2 = – pow(1 – pow((v3+k1/3)/c,2),3/2)*(((v3+k1/3)/mmin)*b + k)*dt;
kr3 = – pow(1 – pow((v3+2*k2/3)/c,2),3/2)*(((v3+2*k2/3)/mmin)*b + k)*dt;
krt = (kr1 + 3*kr3)/4;

// ————————————— RK 5 ——————————————————————-//
k1 = – pow(1 – pow(v5/c,2),3/2)*((v5/mmin)*b + k)*dt;
k2 = – pow(1 – pow((v5+k1)/c,2),3/2)*(((v5+k1)/mmin)*b + k)*dt;
k3 = – pow(1 – pow((v5+(k1+k2)/2)/c,2),3/2)*(((v5+(k1+k2)/2)/mmin)*b + k)*dt;
k4 = – pow(1 – pow((v5+(14*k1+5*k2-3*k3)/64)/c,2),3/2)*(((v5+(14*k1+5*k2-3*k3)/64)/mmin)*b + k)*dt;
k5 = – pow(1 – pow((v5+(-12*k1-12*k2+8*k3+64*k4)/96)/c,2),3/2)*(((v5+(-12*k1-12*k2+8*k3+64*k4)/96)/mmin)*b + k)*dt;
k6 = – pow(1 – pow((v5+(-9*k2 +5*k3 +16*k4 +36*k5)/64)/c,2),3/2)*(((v5+(-9*k2 +5*k3 +16*k4 +36*k5)/64)/mmin)*b + k)*dt;
kt = (7*k1 + 7*k3 +32*k4 +12*k5 +32*k6)/90;

cout << i << “\t” << t << “\t” << mmin << “\t” << v1 << “\t” << v3 << “\t” << v5 << “\t” <<endl;
file << i << “\t” << t << “\t” << mmin << “\t” << v1 << “\t” << v3 << “\t” << v5 << “\t” <<endl;

v1 += k;
v3 += krt;
v5 += kt;

t += dt;
i ++;
}
return 0;
}

================================== Gráfica ==================================

grafica

Advertisements

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s