Runge kuta
Enviado por María José Andaur Estica • 29 de Julio de 2017 • Documentos de Investigación • 6.784 Palabras (28 Páginas) • 281 Visitas
> restart;#
#
# Escribir la función a la cual se le quiere aplicar el método de Runga
# Kutta.
#
> f:=(x,y)->(.313*29.51)*exp(29.51*x)/((1670*x*9.8)*(1-(173.611*(1/.15-1
> /.25))/(950.655*(2670*0.28e-1)*(1-.15)^26.97)));
>
>
>
f := (x, y) ->
0.313 (29.51) exp(29.51 x)
-----------------------------------------------------
/ / 1 1 \ \
| 173.611 |---- - ----| |
| \0.15 0.25/ |
1670 x 9.8 |1 - ------------------------------------|
| 26.97|
\ 950.655 2670 (0.028) (1 - 0.15) /
> x0:=0.15;# Condición inicial (Concentración de entrada en el
# espesador).
x0 := 0.15
> y0:=0; # Altura inicial del espesador.
y0 := 0
> xf:=0.25;# Condición final (Concentración de descarga en el
# espesador).
xf := 0.25
> soln:=(dsolve({ODE,y(x0)=y0})):
Error, (in dsolve) not an ODE system with respect to the unknowns
{y(3/20)}
> h:=(xf-x0)/10; # Intervalos en los que se va aplicar el método o
# número de iteraciónes, en este caso 10.
h := 0.01000000000
> X[0]:=x0: Y[0]:=y0: # Condiciones iniciales.
> X[1]:=X[0]+h;# Calculo de la concentración
# correspondiente al primer intervalo o primera iteración.
#
X[1] := 0.1600000000
> k1:=f(X[0],Y[0]);# Calculo de las K d runge Kuta para la
# primera iteración.
>
> k2:=f(X[0]+h,Y[0]+0.5*k1*h);
> k3:=f(X[0]+h,Y[0]+0.67*k2*h);
> k4:=f(X[0]+h,Y[0]+k3*h);
> Y[1]:=Y[0]+0.125*(k1+3*k2+3*k3+k4)*h;
k1 := 0.6579437548
k2 := 0.8285531011
k3 := 0.8285531011
k4 := 0.8285531011
Y[1] := 0.008072269328
> X[0];
> #
> Y[0];
> k1:=f(X[0],Y[0]);
0.15
0
k1 := 0.6579437548
> X[0]+h;
> Y[0]+0.5*k1*h;
> k2:=f(X[0]+h,Y[0]+0.5*k1*h);
0.1600000000
0.003289718774
k2 := 0.8285531011
> X[0]+h;
> Y[0]+0.67*k2*h;
> k3:=f(X[0]+h,Y[0]+0.67*k2*h);
0.1600000000
0.005551305777
k3 := 0.8285531011
> X[0]+h;
> Y[0]+k3*h;
> k4:=f(X[0]+h,Y[0]+k3*h);
0.1600000000
0.008285531011
k4 := 0.8285531011
> 0.125*(k1+3*k2+3*k3+k4);
> Y[1]:=Y[0]+0.125*(k1+3*k2+3*k3+k4)*h; # Altura del espesador en la
# primera iteración.
0.8072269327
Y[1] := 0.008072269328
> XY:=[[X[0],Y[0]],[X[1],Y[1]]]:
>
> X[2]:=X[1]+h;# Calculo de la concentración correspondiente al
# segundo intervalo o segunda iteración.
X[2] := 0.1700000000
> k1:=f(X[1],Y[1]);
> k2:=f(X[1]+h,Y[1]+0.5*k1*h);
> k3:=f(X[1]+h,Y[1]+0.67*k2*h);
> k4:=f(X[1]+h,Y[1]+k3*h);
> Y[2]:=Y[1]+0.125*(k1+3*k2+3*k3+k4)*h;
k1 := 0.8285531011
k2 := 1.047494400
k3 := 1.047494400
k4 := 1.047494400
Y[2] := 0.01827353671
> X[1];
> Y[1];
> k1:=f(X[1],Y[1]);
0.1600000000
0.008072269328
k1 := 0.8285531011
> X[1]+h;
> Y[1]+0.5*k1*h;
> k2:=f(X[1]+h,Y[1]+0.5*k1*h);
0.1700000000
0.01221503483
k2 := 1.047494400
> X[1]+h;
> Y[1]+0.67*k2*h;
> k3:=f(X[1]+h,Y[1]+0.67*k2*h);
0.1700000000
0.01509048181
k3 := 1.047494400
> X[1]+h;
> Y[1]+k3*h;
> k4:=f(X[1]+h,Y[1]+k3*h);
0.1700000000
0.01854721333
k4 := 1.047494400
> 0.125*(k1+3*k2+3*k3+k4);
> Y[2]:=Y[1]+0.125*(k1+3*k2+3*k3+k4)*h;# Altura del espesador
# en la segunda iteración.
1.020126738
Y[2] := 0.01827353671
> XY:=[[X[0],Y[0]],[X[1],Y[1]],[X[2],Y[2]]]:
>
> X[3]:=X[2]+h; # Calculo de la concentración correspondiente al
# tercer intervalo o tercera iteración.
X[3] := 0.1800000000
> k1:=f(X[2],Y[2]);
> k2:=f(X[2]+h,Y[2]+0.5*k1*h);
> k3:=f(X[2]+h,Y[2]+0.67*k2*h);
> k4:=f(X[2]+h,Y[2]+k3*h);
> Y[3]:=Y[2]+0.125*(k1+3*k2+3*k3+k4)*h;
...