Método de Euler Sencillo en Scilab, múltiples ecuaciones.

Método de Euler Sencillo en Scilab

Autor : Marco Polo Jácome Toss	
Fecha:  05 de Enero del 2022
Licencia : GNU General Public License (GPL) 3.0
Plataforma : Scilab
Code : Euler Sencillo.

Para dar solución a ecuaciones diferenciales ordinarias de la forma :

$$ \frac{dy}{dx}=f\left ( x, y \right ) $$

El método de Euler consiste en encontrar iterativamente la solución de una ecuación diferencial de primer orden (EDO) y valores iniciales conocidos para un rango de valores. Partiendo de un valor inicial "Xo" y avanzando con un paso "h", se pueden obtener la solución de la siguiente manera:

$$ Y_{k+1}=Y_{k}+hf\left (x_{k},y_{k} \right ) $$

donde Y es la solución de la EDO y f es la ecuación diferencial en función de las variables independiente.

Scilab Método de Euler Sencillo.

El método de Euler Sencillo puede ser fácilmente programado en Scilab para múltiples ecuaciones diferenciales.

Primero se define la función que contendrá las ecuaciones diferenciales.

Se ha definido una EDO del tipo :

$$ \frac{dy}{dt}=f\left ( t, y \right ) $$

Algo confuso para algunos cuando no se ve la variable equis.

Segundo paso se define la ecuación el Método de Euler sencillo, observe cuidadosamente que se tiene como salida [y, x, t] y se pide (t0, tf, inicond_x, inicond_y). Para poder hacer uso de esta función es necesario definir el tiempo inicial y final de simulación como las condiciones iniciales.

Planteada la ecuación diferencial:

$$ \frac{dy}{dt}=y-t \: \: \: y\left (0 \right ) =2 $$

Se realiza el planteamiento de las condiciones iniciales:

El paso se saca mediante el tiempo inicial y final entre el número de muestras:

$$ h=\frac{t_{f}-t_{i}}{n} $$

Ahora bien, si incrementamos el número de muestras el paso será más pequeño y se obtendrá una solución más optima. Haciendo el llamado de la función:

[y, x, t] = eulersencillo(t0, tf, n, inicond_x, inicond_y)
disp([t,y'])

Graficando la solución con Scilab mediante el bloque siguiente:

Gráfica, solución EDO

Graficando se obtiene lo siguiente:

Solución de ecuación diferencial ordinaria, circuito RC con fuente de CD

Solución de ecuación diferencial con ODE de Scilab

ODE (Ordinary Differential Equation Solver) es una rutina para resolver ecuaciones diferenciales ordinarias de esta manera usted puede dar solución rápida a una Ecuación diferencial que se encuentre despejada debidamente.

En este ejemplo veremos como dar solución a una ecuación diferencial ordinaria cuando se alienta un circuito RL en serie con una fuente de tensión de directa donde se estudiarán dos casos prácticos, uno sencillo donde únicamente se plasma la EDO y el otro será necesario usar una estrategia distinta usando un poco de programación propia de Scilab.

Caso 1 : Circuito RC Serie con fuente de corriente directa, conectar

En la figura 1 se muestra un circuito RC en serie conectado a una batería de 2000V (E), el valor del capacitor ( C ) y resistencia ( R ) son de 2000 Ω y 40 μF respectivamente. Este circuito al cerrar el interruptor la corriente empezará a fluir hasta llegar a cargar el capacitor de esta manera la corriente disminuirá hasta ser cero.

La EDO (Ecuación Diferencial Ordinaria) del circuito es :

$$ \large E=Ri+\frac{q}{C} $$
La integral de la carga esta dada por :
$$ \large q = \int i dt $$
Sustituyendo en la EDO (Ecuación Diferencial Ordinaria) del circuito se obtiene :
$$ \large E=Ri+\frac{1}{C}\int i dt $$
Diferenciando con respecto al tiempo se obtiene :
$$ \large 0=R\frac{di}{dt}+\frac{i}{C} $$
Despejando el diferencial de la corriente :
$$ \large \frac{di}{dt}=- \frac{i}{CR} $$
Ahora bien, necesitamos plasmar las condiciones iniciales que son :

• Al cerrar el interruptor la corriente iniciará su recorrido en la malla, por lo tanto la corriente inicial es E/R.
• El tiempo inicial es cero de esta manera al cerrar el interruptor empezará a correr el tiempo.
• El tiempo final será necesario para finalizar la simulación del circuito, usaremos 0.5 segundos.

Debemos entender que primero se obtiene la corriente eléctrica y luego esta corriente se sustituye en la EDO que no incluye la integral, despejando la carga se obtiene :
$$ \large q=C(E-Ri) $$
El bloque de código siguiente describe todo el proceso para obtener la corriente y la carga, este código es un poco distinto pero se obtiene primero la corriente eléctrica y después la carga.

global C R E

C=0.00004;
R=2000;
E=200;
    
function ydot = EcuDif(t, y)
    
    ydot=-y/(C*R)
    
endfunction

function q =Carga(y)
    
    q=C*(E-(y*R))
    
endfunction


y0=200/2000;
t0=0.0;
t=0:0.01:0.5;
y = ode(y0,t0,t,EcuDif);
q =Carga(y)


 
scf(1)
plot(t,y,'Color','r','LineWidth',3)
legend('I[Amperios]')
xlabel('Tiempo (Segundos)','fontsize',2)
ylabel('Corriente [Amperios]','fontsize',2)
title('Corriente circuito RC con fuente de CD','fontsize',3)
xgrid



scf(2)
plot(t,q,'Color','b','LineWidth',3)
legend('q[Culombios]')
xlabel('Tiempo (Segundos)','fontsize',2)
ylabel('Carga q [Culombios]','fontsize',2)
title('Carga en circuito RC con fuente de CD','fontsize',3)
xgrid

La rutina ODE nos arroja la respuesta de la ecuación diferencial, una respuesta única como se muestra :

y = ode(y0,t0,t,f);

Corriendo la simulación se obtiene la corriente :

La carga del capacitor es :

Si analizamos un poco más el capacitor se ha cargado en un tiempo de 0.5 segundos. Si dividimos el resultado final de la carga obtenida entre la tensión dará como resultado los microfaradios del capacitor.

$$ \LARGE C = \frac{q}{V}=\frac{0.008\: Columbios}{120\: Volts }=40\mu F $$

Solución de ecuación diferencial ordinaria, circuito RL con fuente de CD

Solución con ODE de Scilab

ODE (Ordinary Differential Equation Solver) es una rutina para resolver ecuaciones diferenciales ordinarias de esta manera usted puede dar solución rápida a una Ecuación diferencial que se encuentre despejada debidamente.

En este ejemplo veremos como dar solución a una ecuación diferencial ordinaria cuando se alimenta un circuito RL en serie con una fuente de tensión de directa donde se estudiarán dos casos prácticos, uno sencillo donde únicamente se plasma la EDO y el otro será necesario usar una estrategia distinta usando un poco de programación propia de Scilab.

Caso 1 : Circuito RL Serie con fuente de corriente directa, conectar

En la figura 1 se muestra un circuito RL en serie conectado a una batería de 10V (E), el valor de la resistencia ® e inductancia (L) son de 20 Ω y 0.6H respectivamente. Este circuito al cerrar el interruptor la corriente empezará a fluir y se incrementará hasta llegar a una corriente constante a lo largo del tiempo por este motivo se establece un tiempo final de simulación.

La EDO (Ecuación Diferencial Ordinaria) del circuito es :

$$ \large E=Ri+L\frac{di}{dt} $$
Dejando libre la derivada de la corriente con respecto al tiempo se obtiene :
$$ \large \frac{di}{dt}=\frac{E-Ri}{L} $$
Ahora bien, necesitamos plasmar las condiciones iniciales que son :

• Al cerrar el interruptor la corriente iniciará su recorrido en la malla, por lo tanto la corriente inicial es cero.
• El tiempo inicial es cero de esta manera al cerrar el interruptor empezará a correr el tiempo.
• El tiempo final será necesario para finalizar la simulación del circuito, usaremos 0.5 segundos.

El siguiente código describe la función que contiene la ecuación diferencial ordinaria y las condiciones iniciales.

function ydot=EcuDif(t, y)
    L=0.6;
    R=20;
    E=10;
    ydot=(E-R*y)/L 
endfunction

t0=0.0;
t=0:0.001:0.5;
y = ode(y0,t0,t,EcuDif); 
scf(1)
plot(t,y','-')
legend('I[A]')
xlabel('Tiempo (Segundos)','fontsize',2)
ylabel('Corriente [Amperios]','fontsize',2)
title('Corriente circuito RL con fuente de CD','fontsize',2)
xgrid

Finalmente se obtiene la figura 1 que se muestra a continuación :

Caso 2: Circuito RL Serie con fuente de corriente directa, desconectar

En la figura 2 se muestra un circuito RL en serie conectado a una batería de 10V (E), el valor de la resistencia ® e inductancia (L) son de 20 Ω y 0.6H respectivamente. Este circuito en particular ha alcanzado su corriente nominal pero el interruptor se abre de esta manera la corriente empezará a disminuir drásticamente hasta ser igual a cero.

La EDO (Ecuación Diferencial Ordinaria) del circuito es :

$$ \large E=Ri+L\frac{di}{dt} $$
Dejando libre la derivada de la corriente con respecto al tiempo se obtiene :
$$ \large \frac{di}{dt}=\frac{E-Ri}{L} $$
Ahora bien, necesitamos plasmar las condiciones iniciales que son :

• Al abrir el interruptor la corriente descenderá en la malla, por lo tanto la corriente al final será cero.
• El tiempo inicial es cero de esta manera al abrir el interruptor empezará a correr el tiempo.
• El momento cuando se abrirá el interruptor, usaremos 0.2 segundos.
• El tiempo final será necesario para finalizar la simulación del circuito, usaremos 0.5 segundos.

El siguiente código describe la función que contiene la ecuación diferencial ordinaria y las condiciones iniciales.

function ydot=EcuDiff(t, y)
    L=0.6;
    R=20;
    E=10;
    if t<0.2 then
        E=10;
        ydot=(E-R*y)/L       
    else
        E=0;
        ydot=(E-R*y)/L 
    end
    
endfunction

t0=0.0;
t=0:0.001:0.5;

for i=1:max(size(t))
    
    if t(i)<0.2 then
        y0=0.5
        y(i) = ode(y0,t0,t(i),EcuDiff);        
    else
        y0=0
        y(i) = ode(y0,t0,t(i),EcuDiff); 
    end
    
end   
scf(1)
plot(t,y','-')
legend('I[A]')
xlabel('Tiempo (Segundos)','fontsize',2)
ylabel('Corriente [Amperios]','fontsize',2)
title('Corriente circuito RL con fuente de CD','fontsize',2)
xgrid

Finalmente se obtiene la figura 2 que se muestra a continuación :

Explicación del Caso 2

La simulación de este segundo caso es un poco complicada debido a que necesitamos programar ciertas condiciones. La función de la ecuación diferencial debe poseer una condición importante, la desconexión del interruptor, entonces usamos la condición if.

La condición es sencilla cuando se conoce el tiempo en le momento que abre el interruptor en nuestro caso 0.2 segundos. Después de ese tiempo la tensión será cero por lo tanto se coloca de nuevo la ecuación diferencial con esta nueva tensión.

function ydot=EcuDiff(t, y)
    L=0.6;
    R=20;
    E=10;
    //Tiempo de apertura
    if t<0.2 then
        E=10;
        ydot=(E-R*y)/L 
    //Depués de la apertura del interruptor    
    else
        E=0;
        ydot=(E-R*y)/L 
    end
    
endfunction

Para poder correr esta condición usaremos el ciclo for de esta manera se plantea la misma condición, lo importante aquí es analizar la cantidad de muestras por lo tanto necesitamos del ciclo for, finalmente las variables van cambiando conforme corre el tiempo siendo necesario actualizarlas y para realizar esto en Scilab usamos y(i) donde i es un número del total de muestras. La misma situación se plantea con el tiempo debido a que va cambiando y te preguntarás por qué no usamos t(i) ni y(i) en la función donde definimos la ecuación diferencial; la respuesta es muy sencilla, actualizamos el tiempo fuera de la función y no dentro haciendo el trabajo de la simulación mucho más sencillo.

for i=1:max(size(t))
    
    if t(i)<0.2 then
        y0=0.5
        y(i) = ode(y0,t0,t(i),EcuDiff);        
    else
        y0=0
        y(i) = ode(y0,t0,t(i),EcuDiff); 
    end
    
end 

Métodos numéricos, Método de la Bisección para obtener raíces con Scilab

Método de la bisección

El método de bisección, también conocido como corte binario, de participación de intervalos o de Bolzano, es una técnica de búsqueda incremental que consiste en dividir un intervalo sucesivamente a la mitad.

Su funcionamiento se basa en la evaluación del valor de la función en el punto medio del intervalo. Si la función cambia de signo en dicho punto, se determina que la raíz se encuentra en uno de los subintervalos resultantes. El proceso se repite en el subintervalo que contiene el cambio de signo, dividiéndolo nuevamente a la mitad y evaluando la función en el nuevo punto medio.

Este proceso iterativo se continúa hasta que se alcanza la precisión deseada en la aproximación a la raíz.

Las principales ventajas del método de bisección son:

• Simpleza y facilidad de comprensión.
• No requiere información adicional sobre la función a evaluar.
• Convergencia garantizada en la mayoría de los casos.

Sin embargo, también presenta algunas desventajas:

• Convergencia lenta en algunos casos.
• Precisión limitada en comparación con otros métodos numéricos.
• No es aplicable para la búsqueda de raíces complejas.

En este ejemplo vamos a obtener las posibles raíces de la función siguiente :

Para poder resolverlo se usará el bloque de código siguiente :

Autor : Marco Polo Jácome Toss  
Fecha de creación :  29 de Agosto del 2021
Licencia : GNU General Public License (GPL) 3.0
Plataforma : Scilab, Método de la Bisección

La siguiente función es el método de la bisección.

// fun = function handle
// x = Intervalo de una posible raíz
// tol = Máxima tolerancia
// maxit = Máximo número de iteraciones
// Root: Raíz

clear; 
format("v",8)
clc
clear

function [root] = Bisection(fun,x,tol,maxit)

if fun(x(1)) > 0 then
    xu = x(1);    xl = x(2);
else
    xu = x(2);    xl = x(1);
end

Ea = 1;
iter = 1;

while(1)
    
    xr(iter) = (xl(iter) + xu(iter)) / 2;

    if fun(xr(iter)) > 0 then
        xu(iter+1) = xr(iter);
        xl(iter+1) = xl(iter);
    elseif fun(xr(iter)) < 0 then
        xl(iter+1) = xr(iter);
        xu(iter+1) = xu(iter);
    else
        break
    end
    
    if iter>1 then
        Ea(iter) = 100 * abs((xr(iter) - xr(iter-1)) / xr(iter));
    end
    
    disp([iter xr(iter) Ea(iter) xu(iter+1) xl(iter)])

    if Ea(iter) < tol | iter == maxit then
        break
    end
    iter = iter + 1;
end
root = xr(iter);
endfunction

La función de este ejemplo debe escribirse de la forma siguiente :

function f = fun1(x)
    f = 3*x-(x*x)-(5*x*x*x*x*x)
endfunction

Las variables iniciales para la función son :

x = [-5 0];
tol = 1e-3;
maxit = 100;
root = Bisection(fun1,x,tol,maxit);
disp("Raíz = ", root)

//Segundo intervalo

x = [0.5 1];
tol = 1e-3;
maxit = 100;
root = Bisection(fun1,x,tol,maxit);
disp("Raíz = ", root)

Graficando la función

Graficando la función se puede observar dos raíces entre el intervalo -1 a 1.

Raíces de la función

La raíz en el intervalo de -5 a 0 es :

    1.  - 2.5    1.  
    2.  - 1.25    100.  
    3.  - 0.625    100.  
    4.  - 0.9375    33.333333  
    5.  - 1.09375    14.285714  
    6.  - 1.015625    7.6923077  
    7.  - 0.9765625    4.  
    8.  - 0.9570313    2.0408163  
    9.  - 0.9472656    1.0309278  
    10.  - 0.9423828    0.5181347  
    11.  - 0.9399414    0.2597403  
    12.  - 0.9411621    0.1297017  
    13.  - 0.9417725    0.0648088  
    14.  - 0.9420776    0.0323939  
    15.  - 0.9422302    0.0161943  
    16.  - 0.9423065    0.0080965  
    17.  - 0.9423447    0.0040481  
    18.  - 0.9423256    0.0020241  
    19.  - 0.9423161    0.0010121  
    20.  - 0.9423113    0.0005060  
 root =    
  - 0.9423113

La posible raíz en el intervalo 0 a 1 es :

   1.   0.55   1.   0.55   1.
   2.   0.775   29.0323   0.775   1.
   3.   0.8875   12.6761   0.775   1.
   4.   0.83125   6.76692   0.775   0.8875
   5.   0.80313   3.50195   0.80313   0.83125
   6.   0.81719   1.72084   0.80313   0.83125
   7.   0.81016   0.86789   0.81016   0.81719
   8.   0.81367   0.43207   0.81016   0.81719
   9.   0.81191   0.2165   0.81191   0.81367
   10.   0.81279   0.10813   0.81279   0.81367
   11.   0.81323   0.05404   0.81279   0.81367
   12.   0.81301   0.02703   0.81301   0.81323
   13.   0.81312   0.01351   0.81312   0.81323
   14.   0.81318   0.00676   0.81318   0.81323
   15.   0.8132   0.00338   0.8132   0.81323
   16.   0.81322   0.00169   0.81322   0.81323
   17.   0.81323   0.00084   0.81322   0.81323

--> disp("Raz = ", root)

  "Raíz = "

   0.81323

Métodos numéricos, Método de la Bisección para obtener raíces con Scilab por Marco Polo Jácome Toss se distribuye bajo una Licencia Creative Commons Atribución-CompartirIgual 4.0 Internacional.

Archivos CSV con Scilab

Scilab : Guardar resultados en un archivo CSV

En este pequeño tutorial les explicaré como guardar los resultados obtenidos en Scilab en un archivo de Excel.

Primero generamos algunos resultados mediante el código siguiente :

fhz=60;
ciclos=2;
tfinal=ciclos/fhz;
t=0:0.0001:tfinal;
wt=2*%pi()*fhz*t
y=sin(wt)
scf(1)
plot(t,y,'r')
title("Función senoidal")
xlabel("X")
ylabel("Y")
legend("y=sin(377t)")
xgrid(35)

El código anterior genera la gráfica :

La gráfica anterior puede ser editada en el menú editar en propiedades de la figura.

Segundo capturamos la ruta donde se encuentra el archivo *.sce mediante

//pwd() sirve para obtener la ruta donde se encuentra el archivo SCilab que ejecutamos.
ruta=strcat([pwd(),'\datos_obtenidos.csv'])

Las variables se acomodan mediante un arreglo matricial para ser guardado en el archivo de Excel.

--> [t',y']
 ans  =

   0.       0.       
   0.0001   0.0376902
   0.0002   0.0753268
   0.0003   0.1128564
   0.0004   0.1502256
   0.0005   0.1873813
   0.0006   0.2242708
   0.0007   0.2608415
   0.0008   0.2970416
   ......   .........

Generamos el archivo CSV en la ruta donde fue ejecutado el archivo de Scilab.

//*Los resultados obtenidos se acomodan en columnas y escribimos el archivo *.CSV
write_csv([t',y'],ruta)

Abrimos el archivo CSV

Finalmente el código completo

Guarde las líneas de código con extensión *.sce en la carpeta Documentos o en la que guste y corra el archivo para obtener los resultados obtenidos en esta publicación.

fhz=60;
ciclos=2;
tfinal=ciclos/fhz;
t=0:0.0001:tfinal;
wt=2*%pi()*fhz*t
y=sin(wt)
scf(1)
plot(t,y,'r')
title("Función senoidal")
xlabel("X")
ylabel("Y")
legend("y=sin(377t)")
xgrid(35)

ruta=strcat([pwd(),'\datos_obtenidos.csv'])
write_csv([t',y'],ruta)

Métodos de integración de ecuaciones diferenciales

Métodos de integración Runge-Kutta y Regla Trapezoidal

El método de Runge-Kutta (RK) es un conjunto de métodos iterativos (implícitos y explícitos) para la aproximación de soluciones de ecuaciones diferenciales ordinarias, en esta documentación se muestra el de cuarto y sexto orden. Otro método agregado es la Regla Trapezoidal que puede ser explicita, para diferenciar los métodos de solución explicito e implícito se agregan los dos puntos siguientes:

Autor : Marco Polo Jácome Toss	
Fecha de creación :  12 de Septiembre del 2017
Licencia : GNU General Public License (GPL) 3.0
Plataforma : Scilab
Código fuente creado para dar solución numérica a múltiples ecuaciones diferenciales.
Actualizado : 04 de Abril del 2022

Método de integración explicito. En este método es posible calcular la aproximación en cada paso directamente evaluando la función f(x,y), como ejemplo de un método de integración explicito se muestra la regla del punto medio.

Donde

Método de integración implicito. Las aproximaciones en este método vienen definidas por un sistema de ecuaciones implícito. En la siguiente ecuación se muestra la forma implícita del método de integración del punto medio.

La ventaja del "scripting" es poder dar solución en forma ordenada a n número de ecuaciones diferenciales.

Problemas de valor inicial

Se considera el problema de valor inicial (PVI) de ecuaciones diferenciales ordinarias (EDO) :

Lo anterior es bastante importante debido a que debemos indicar las condiciones iniciales para poder dar solución a una EDO.

1.1 ¿Por qué utilizar este código fuente ?

Debido a su simplicidad de resolver varias ecuaciones de estado es flexible para dar solución a modelos de máquinas eléctricas como para otros tipos de modelos que dependan principalmente del tiempo o en caso contrario que no involucre esta variable.

Las líneas de código siguientes pertenecen al archivo start.sce el cual contiene tres métodos de solución Rk4,Rk3 y Trapezoidal ya antes mencionados.

1.2 Lista de archivos dependientes

La siguiente lista de archivos (dependientes) muestra como se encuentra estructurado en forma general el programa y siempre debe ejecutar el archivo start.sce para poder llamar la lista restante de archivos con extensión *.sci .

1. Start
   • attributes.sci
   • ecuDif.sci
   • rk4.sci (Multi ecuaciones)
   • rk6.sci (Multi ecuaciones)
   • rtrapezoidal (Multi ecuaciones)
   • mnecudif

1.3 Ejecución del código fuente

De manera sencilla iniciamos el Start.sce e inmediatamente se indica el tipo de método a ocupar.

**Seleccione la forma de solución : RK4(1) / RK6(2) / RTRAPEZOIDAL(3) : **

1.4 Ecuaciones diferenciales : ecuDif.sci

En el archivo ecuDif es necesario especificar las ecuaciones de estado, por lo que Xdot contiene dos ecuaciones de estado por este motivo se tiene dos variables Xdot(1,1) y Xdot(2,1). Estas variables se pasan a los archivos rk4.sci y rk6.sci junto con la variable de tiempo para poder dar solución a la EDO.

En este archivo se introducen las constantes y valores necesarios de las ecuaciones de estado con sus respectivas variables de estado representados con Xdot, todo es un arreglo matricial.

1.5 Método de Runge-Kutta 6to Orden : rk6.sci

Es necesario cambiar la línea del archivo rk4 o rk6 que contiene disp([t(i),r(1,i)',r(2,i)']) al incrementar las variables y ecuaciones o únicamente mostrar una al igual que las variables iniciales mostradas en el archivo attributes.sci.

1.6 Método de Runge Kuta 4to Orden : rk4.sci

La función rk4 se muestra en el bloque de código siguiente:

1.7 Método Trapezoidal : rtrapezoidal.sci

1.8 Valores Iniciales de Simulación : attributes.sci

El archivo contiene los atributos de la simulación y debe establecer el número de condiciones iniciales, las cuales dependen de las variables de estado a resolver siendo necesario modificar la impresión disp([t(i),r(1,i)',r(2,i)']) de los archivos rk4 y rk6.

• ti: tiempo inicial de simulación.
• tf: tiempo final de simulación.
• varIni: variables iniciales.
• muestras: muestras.

Los valores iniciales mostrados en el bloque siguiente son el tiempo inicial, tiempo final, valores iniciales en este caso para dos variables y el número de muestras.

global ti tf varIni muestras    
ti=0;
tf=10;
varIni=[0.01, 0.02]; //Cambiar #Xdot
muestras=10000;

1.9 Opciones de solución: mnecudif.sci

Las opciones mencionadas en la parte 1.1 pertenecen al código siguiente :

Ejemplo #1 Ecuación del péndulo

Un péndulo simple se define como una partícula de masa "m" suspendida "O" por un hilo inextensible de longitud "L" y de masa despreciable.

Aplicando la segunda Ley de Newton, siendo la aceleración de la partícula hacia adentro de la trayectoria circular se obtiene:

Siendo la aceleración de la partícula hacia adentro dada por:

La aceleración tangencial de la partícula es:

De la segunda Ley de Newton incluyendo la fricción se obtiene:

El valor de k (fricción) esta dado por :

Las variables de estado son :

Por lo tanto, las ecuaciones de estado son las siguientes:

Estableciendo la ecuaciones en el archivo ecuDif.sci y valores correspondientes.

function [Xdot]=ecuDif(t,x)
  m=0.5;
  b=0.1;
  L=1.5;
  g=9.81; 
  k=b/(m*L);
 Xdot=zeros(2,1); //#Xdot
 Xdot(1,1)=x(2);
 Xdot(2,1)=-(g/L)*sin(x(1))-(k/m)*x(2);
endfunction

La solución paraXdot(1) y Xdot(2) de las ecuaciones de estado es:

Ejemplo #2 Solución analítica y numérica

Utilizando el método de la transformada de Laplace se resolverá la ecuación diferencial siguiente:

 

Aplicando la transformada de Laplace a la Ecuación diferencial anterior se obtiene la solución para Y(S).

   

La solución en el dominio del tiempo se obtiene con la transformada inversa.

   

El resultado anterior se grafica con Scilab mediante el bloque siguiente:

t=0:0.01:1;
y=8*exp(-3*t)-2*cos(2*t)+3*sin(2*t)
plot(t,y)
xlabel("Tiempo, seg.")
ylabel('y(t)')
title('Solución de la Ecuación Diferencial')
xgrid()

La solución de la ecuación diferencial y gráfica se muestra a continuación:

 

Para dar solución es necesario configurar el archivo ecuDif de la forma siguiente:

function [Xdot]=ecuDif(t,x)
 Xdot=zeros(1,1); //#Xdot
 Xdot(1,1)=13*sin(2*t)-3*x(1)
endfunction

Modificamos las condiciones iniciales del archivo attributes como se muestra:

global ti tf varIni muestras    
ti=0;
tf=1;
varIni=[6,0]; //Cambiar #Xdot
muestras=1000;

Se debe ajustar el archivo mnecudif para mostrar únicamente una solución, esta ecuación diferencial tiene una única variable.

function mnecudif(alg)
if(alg == 1) then
[t,r]=rk4(ti,tf,muestras,varIni);
plot(t,[r(1,:)])
legend('Xdot(1)')
xlabel('tiempo, seg')
ylabel('Y')
title('Solución Runge Kutta orden 4')
xgrid
elseif ( alg == 2) then
[t,r]=rk6(ti,tf,muestras,varIni);
plot(t,[r(1,:)'])
legend('Xdot(1)')
xlabel('tiempo, seg')
ylabel('Ysol')
title('Solución Runge Kutta orden 6')
xgrid
elseif (alg == 3) then
[t,r]=rtrapezoidal(ti,tf,muestras,varIni);
plot(t,[r(1,:)'])
legend('Xdot(1)')
xlabel('tiempo, seg')
ylabel('Y')
title('Solución Regla Trapezoidal Explicita')
xgrid
else
disp('Error: Seleccione un método de solución ');
end
endfunction 

Finalmente al realizar estas modificaciones podrás observar la solución siguiente:

Copyright

Copyright © 2017 en adelante, Marco Polo Jácome Toss (https://jacometoss.github.io/Runge_Kutta_4th_6th/). Este programa es software libre: usted puede redistribuirlo y /o modificarlo bajo los términos de la Licencia General GNU (GNU General Public License) publicado por la Fundación para el Software Libre para la versión 3 de dicha Licencia o anterior, o cualquier versión posterior.

Este programa se distribuye con la esperanza de que sea útil pero sin ninguna garantía; incluso sin la garantía implícita de comercialización o idoneidad para un propósito en particular.

Vea la información de Licencia de RK4 para más detalle.

Información de registro