11/Mayo/2015

Instituto Politcnico
Nacional

Alumno:

Ramos Albarrn Fernando

Profesor: Dr. David Romero Romero.

Curso: Anlisis de Sistemas Lineales.

Seccin de Estudios de Posgrado e Investigacin

Escuela Superior de Ingeniera Mecnica y Elctrica.

Sistemas Lineales.

[PROYECTO: SOLUCIN NUMRICA

DE ECUACIONES DIFERENCIALES]

Seccin de Estudios de Posgrado e Investigacin

1

Escuela Superior de Ingeniera Mecnica y Elctrica.

Sistemas Lineales.

Tabla de contenido
INTRODUCCIN................................................................................................ 2
Mtodos Explcitos........................................................................................... 3
Mtodo Predictor - Corrector............................................................................ 3
Mtodos Implcitos........................................................................................... 4
DESARROLLO.................................................................................................... 4
EJERCICIO 1...................................................................................................... 4
A) SOLUCIN POR EL MTODO DE EULER...................................................4
B) SOLUCIN POR EL MTODO TRAPEZOIDAL.............................................6
C) SOLUCIN POR EL MTODO DE RUNGE KUTTA DE 4TO ORDEN..............8
D) SOLUCIN POR EL MTODO DE GEAR DE 1ER ORDEN O EULER
REGRESIVO................................................................................................... 9

EJERCICIO 2.................................................................................................... 11
A) SOLUCIN POR EL MTODO DE HEUN..................................................11
B) SOLUCIN POR EL MTODO TRAPEZOIDAL...........................................13
EJERCICIO 3................................................................................................. 15
A) SOLUCIN POR EL MTODO DE EULER.................................................15
B) SOLUCIN POR EL MTODO TRAPEZOIDAL...........................................17
C) SOLUCIN POR EL MTODO DE RUNGE KUTTA DE 4TO ORDEN............19
D) SOLUCIN POR EL MTODO DE GEAR DE 1ER ORDEN O EULER
REGRESIVO................................................................................................. 20
EJERCICIO 4.................................................................................................... 22
A) DETERMINAR h MXIMO CON EL MTODO DE GEAR DE 1ER ORDEN. . .22
B Y C) SOLUCIN CON EL MTODO DE EULER PARA 2hmax y hmax/2........23
CDIGO DE LOS MTODOS EN FORTRAN 90..............................................25
CONLCUSIONES.............................................................................................. 33
BIBLIOGRAFA................................................................................................. 33

INTRODUCCIN
Sea una ecuacin diferencial ordinaria de la forma
Seccin de Estudios de Posgrado e Investigacin
2

Escuela Superior de Ingeniera Mecnica y Elctrica.

Sistemas Lineales.

(1)
Y sea

(2)
La ecuacin de pendiente ordenada al origen. De acuerdo con esta ecuacin, la
pendiente estimada f se usa para extrapolar desde un valor anterior yi a un
nuevo valor yi+1 en una distancia h (figura 1).
Nuevo valor = valor anterior + pendiente tamao de paso
O bien:
Nuevo valor = valor anterior + funcin evaluada en el valor anterior tamao
de paso

Figura 1 Ilustracin grfica del mtodo de un paso.

Esta frmula se aplica paso a paso para calcular un valor posterior y, por lo
tanto, para trazar la trayectoria de la solucin. Todos los mtodos de un paso
que se expresen de esta forma general, son llamados mtodos explcitos y
tan slo van a diferir en la manera en la que se estima la pendiente. En otras
palabras, se toma la pendiente al inicio del intervalo como una aproximacin
de la pendiente promedio sobre todo el intervalo. Tal procedimiento, llamado
mtodo de Euler. Despus se revisa otro mtodo de un paso que emplean otras
formas de estimar la pendiente que dan como resultado predicciones ms
exactas. Todas estas tcnicas en general se conocen como mtodos de RungeKutta. El mtodo de Euler puede ser considerado tambin como el mtodo de
Adams-Bashford de 1er orden.
Muchas veces se tienen sistemas que involucran constantes de tiempo desde
pequeos valores hasta grandes valores, stos son llamados sistemas rgidos, y
para encontrar soluciones para dichos sistemas, es necesario mezclar los
Seccin de Estudios de Posgrado e Investigacin
3

Escuela Superior de Ingeniera Mecnica y Elctrica.

Sistemas Lineales.
mtodos explcitos con los que se conocen como mtodos implcitos, se forma
entonces un mtodo predictor-corrector, se puede usar solamente algn
mtodo implcito para mejorar la aproximacin, la mayor parte del tiempo
stos ltimos conducen a un sistema de ecuaciones no lineales, por lo cual es
un poco ms complicado aplicar el mtodo. La forma de los mtodos implcitos
es de la forma:
Nuevo valor = valor anterior + funcin evaluada en diferentes valores
incluyendo el actual tamao de paso

Los mtodos usados en ste documento son los siguientes:

Mtodos Explcitos
Mtodo de Euler - Adams-Bashford de 1er Ordeny i+1= yi+ (xi , yi )h

(3)

Mtodos de Runge-Kutta de Cuarto Orden

El ms popular de los mtodos RK es el de cuarto orden. Como en el caso de
los procedimientos de segundo orden, hay un nmero infinito de versiones. La
siguiente, es la forma comnmente usada y, por lo tanto, le llamamos mtodo
clsico RK de cuarto orden:

(4 )

Donde

(5)

Mtodo Predictor - Corrector

Mtodos de Heun

(6)

Seccin de Estudios de Posgrado e Investigacin

4

Escuela Superior de Ingeniera Mecnica y Elctrica.

Sistemas Lineales.

La primer frmula de (6) es la predictora (Euler) que corresponde al

mtodo explcito y la segunda ecuacin es el mtodo corrector que
corresponde al mtodo implcito.
Mtodos Implcitos
Regla Trapezoidal Adams Moulton 2do Ordeny i+1= yi+

h
{ ( xi , yi ) +( x i+1 , y i +1)}
2

(7)

Mtodo de Gear de 1er Orden o Euler Regresivo

y i+1= yi+ h( x i +1 , y i+1 )

(8)

DESARROLLO
x =x 2 ( t )

EJERCICIO 1
x ( t 0 )=1

d e t=0 a t=1 para h=0.1, 0.01, 0.001 y 0.0001

A) SOLUCIN POR EL MTODO DE EULER

Cdigo usado para graficar en MatLab, donde t1 a t4 y x1 a x4 son las salidas
del programa en ForTran, para h=0.1, 0.01, 0.001 y 0.0001 respectivamente. El
cdigo en ForTran del programa para el mtodo de Euler est al final de este
documento. Si se requiere comparar otro valor de h se modifica la funcin en
MatLab, aadiendo una x5 y t5 con sus etiquetas correspondientes.
function Euler_vs_analitica(t1,x1,t2,x2,t3,x3,t4,x4)
%Grfica de la solucin numrica h=0.1:
plot(t1,x1,'-r')
hold on
%Grfica de la solucin numrica h=0.01:
plot(t2,x2,'-m')
hold on
%Grfica de la solucin numrica h=0.001:
plot(t3,x3,'-g')
hold on
%Grfica de la solucin numrica h=0.0001:
plot(t4,x4,'-y')
hold on
%Solucin analtica:

Seccin de Estudios de Posgrado e Investigacin

5

Escuela Superior de Ingeniera Mecnica y Elctrica.

Sistemas Lineales.
Q=length(t4);
SolAn=zeros(Q,1);
for i=1:Q
SolAn(i)=-1./(t4(i)-1);
end
plot(t4,SolAn,'-b')
title('Comparacin del Mtodo de Euler con la solucin analtica');
xlabel('Tiempo t (seg)');
ylabel('x(t)');
legend('h=0.1','h=0.01','h=0.001','h=0.0001','Solucin verdadera');
end
Comparacin del Mtodo de Euler con la solucin analtica

3000

h=0.1
h=0.01
h=0.001
h=0.0001
Solucin verdadera

2500

x(t)

2000

1500

1000

500

0.1

0.2

0.3

0.4

0.5
Tiempo t (seg)

0.6

0.7

0.8

0.9

Haciendo zoom en la grfica se observa que para un paso pequeo la

aproximacin es ms exacta para este sistema.
Comparacin del Mtodo de Euler con la solucin analtica
180
h=0.1
h=0.01
h=0.001
h=0.0001
Solucin verdadera

160

140

x(t)

120

100

80

60

40

20

0.94

0.95

0.96

0.97
Tiempo t (seg)

Seccin de Estudios de Posgrado e Investigacin

6

0.98

0.99

Escuela Superior de Ingeniera Mecnica y Elctrica.

Sistemas Lineales.

B) SOLUCIN POR EL MTODO TRAPEZOIDAL

Se usa exactamente el mismo cdigo de MatLab para graficar
mtodo de Euler mostrado anteriormente, donde t1 a t4 y x1 a
salidas del programa en ForTran, para h=0.1, 0.01, 0.001
respectivamente. El cdigo en ForTran del programa para el mtodo
est al final de este documento.

que en el
x4 son las
y 0.0001
trapezoidal

Comparacin del Mtodo Trapezoidal con la solucin analtica

14000

h=0.1
h=0.01
h=0.001
h=0.0001
Solucin verdadera

12000
10000

x(t)

8000
6000
4000
2000
0

0.1

0.2

0.3

0.4
0.5
0.6
Tiempo t (seg)

0.7

0.8

0.9

Haciendo zoom en varias partes de la grfica se puede observar que el mtodo

de la regla trapezoidal es mucho ms exacto que el mtodo de Euler, para
h=0.0001 es casi exacto para ste sistema, incluso para h=0.001 es ms
exacto que el mtodo de Euler.

Seccin de Estudios de Posgrado e Investigacin

7

Escuela Superior de Ingeniera Mecnica y Elctrica.

Sistemas Lineales.
Comparacin del Mtodo Trapezoidal con la solucin analtica

800

h=0.1
h=0.01
h=0.001
h=0.0001
Solucin verdadera

700

600

x(t)

500

400

300

200

100

0.94

0.95

0.96

0.97
Tiempo t (seg)

0.98

0.99

Comparacin del Mtodo Trapezoidal con la solucin analtica

588

h=0.1
h=0.01
h=0.001
h=0.0001
Solucin verdadera

586

x(t)

584

582

580

578

576
0.9976

0.9977

0.9978

0.9979

0.998
0.9981
Tiempo t (seg)

Seccin de Estudios de Posgrado e Investigacin

8

0.9982

0.9983

0.9984

0.9985

Escuela Superior de Ingeniera Mecnica y Elctrica.

Sistemas Lineales.
Comparacin del Mtodo Trapezoidal con la solucin analtica
2400

h=0.1
h=0.01
h=0.001
h=0.0001
Solucin verdadera

2395

2390

x(t)

2385

2380

2375

2370

2365

2360

0.9995

0.9995

0.9996

0.9996
Tiempo t (seg)

0.9996

0.9996

C) SOLUCIN POR EL MTODO DE RUNGE KUTTA DE 4TO

ORDEN
Usando la misma funcin para graficar en MatLab se obtiene la siguiente
grfica.

Seccin de Estudios de Posgrado e Investigacin

9

Escuela Superior de Ingeniera Mecnica y Elctrica.

Sistemas Lineales.
4

Comparacin del Mtodo Runge Kutta de 4to orden con la solucin analtica

x 10

h=0.1
h=0.01
h=0.001
h=0.0001
Solucin verdadera

x(t)

0.1

0.2

0.3

0.4

0.5
Tiempo t (seg)

0.6

0.7

0.8

0.9

Al hacer zoom se observa que, como el sistema no es rgido; se obtiene una

solucin parecida al mtodo trapezoidal y adems; debido a que su pendiente
se fragmenta en ms pedazos ms pequeos que el mtodo de Euler es ms
exacto que ste ltimo.
Comparacin del Mtodo Runge Kutta de 4to orden con la solucin analtica
h=0.1
h=0.01
h=0.001
h=0.0001
Solucin verdadera

1800

1600

1400

1200

x(t)

1000

800

600

400

200

0
0.992

0.993

0.994

0.995

0.996
Tiempo t (seg)

Seccin de Estudios de Posgrado e Investigacin

10

0.997

0.998

0.999

Escuela Superior de Ingeniera Mecnica y Elctrica.

Sistemas Lineales.
Comparacin del Mtodo Runge Kutta de 4to orden con la solucin analtica
984.7

h=0.1
h=0.01
h=0.001
h=0.0001
Solucin verdadera

984.6

984.5

x(t)

984.4

984.3

984.2

984.1

984
0.999

0.999

0.999

0.999

0.999

0.999
Tiempo t (seg)

0.999

0.999

0.999

0.999

D) SOLUCIN POR EL MTODO DE GEAR DE 1ER ORDEN O

EULER REGRESIVO
Usando la misma funcin para graficar en MatLab se obtiene la siguiente
figura.
4

1.5

Comparacin del Mtodo de Gear de 1er orden con la solucin analtica

x 10

h=0.1
h=0.01
h=0.001
h=0.0001
Solucin verdadera

x(t)

0.5

-0.5

-1

0.1

0.2

0.3

0.4

0.5
Tiempo t (seg)

0.6

Seccin de Estudios de Posgrado e Investigacin

11

0.7

0.8

0.9

Escuela Superior de Ingeniera Mecnica y Elctrica.

Sistemas Lineales.

Haciendo zoom en la figura se observa que para valores cercanos a donde la

funcin es discontinua el mtodo muestra inestabilidad, se observa que es ms
exacto que el Euler explcito pero menos exacto que el trapezoidal o incluso el
Runge Kutta de 4to orden.
Comparacin del Mtodo de Gear de 1er orden con la solucin analtica

4000

3000

h=0.1
h=0.01
h=0.001
h=0.0001
Solucin verdadera

x(t)

2000

1000

-1000

0.8

0.85

Tiempo t (seg)

0.9

Seccin de Estudios de Posgrado e Investigacin

12

0.95

Escuela Superior de Ingeniera Mecnica y Elctrica.

Sistemas Lineales.

EJERCICIO 2
3

x =( xt ) +3 t

x ( 0 )=0

d e t=0 a t =1 para h=22 , 24 ,26 , 28 ,210 y 212

A) SOLUCIN POR EL MTODO DE HEUN
Cdigo usado para graficar en MatLab, donde t1 a t6 y x1 a x6 son las salidas
del programa en ForTran, para

10

h=2 , 2 , 2 ,2 , 2

12

y2

respectivamente.

El cdigo en ForTran del programa para el mtodo de Heun est al final de este
documento. Si se requiere comparar otro valor de h se modifica la funcin en
MatLab, aadiendo una x7 y t7 con sus etiquetas correspondientes.
function heun_vs_analitica2(t1,x1,t2,x2,t3,x3,t4,x4,t5,x5,t6,x6)
%Grfica de la solucin numrica h=2^-2:
plot(t1,x1,'-y')
hold on
%Grfica de la solucin numrica h=2^-4:
plot(t2,x2,'-m')
hold on
%Grfica de la solucin numrica h=2^-6:
plot(t3,x3,'-c')
hold on
%Grfica de la solucin numrica h=2^-8:
plot(t4,x4,'-r')
hold on
%Grfica de la solucin numrica h=2^-10:
plot(t5,x5,'-g')
hold on
%Grfica de la solucin numrica h=2^-12:
plot(t6,x6,'-k')
hold on
%Solucin analtica:
Q=length(t6);
SolAn=zeros(Q,1);
for i=1:Q
SolAn(i)=(t6(i)).^3;
end
plot(t6,SolAn,'-b')
title('Comparacin del Mtodo de Heun con la solucin analtica');
xlabel('Tiempo t (seg)');
ylabel('x(t)');
legend('h=2^-2','h=2^-4','h=2^-6','h=2^-8','h=2^-10','h=2^-12','Solucin
verdadera');
end

Seccin de Estudios de Posgrado e Investigacin

13

Escuela Superior de Ingeniera Mecnica y Elctrica.

Sistemas Lineales.
Comparacin del Mtodo de Heun con la solucin analtica

h=2-2

0.9

h=2-4

0.8

h=2-6
h=2-8

0.7

h=2-10

x(t)

0.6

h=2-12
Solucin verdadera

0.5
0.4
0.3
0.2
0.1
0

0.1

0.2

0.3

0.4
0.5
0.6
Tiempo t (seg)

0.7

0.8

0.9

Haciendo zoom en la grfica se observa que para un paso pequeo la

aproximacin es ms exacta para este sistema.
Comparacin del Mtodo de Heun con la solucin analtica

0.32

h=2-2

0.3

h=2-4

0.28

h=2-6

0.26

h=2-8
h=2-10

x(t)

0.24

h=2-12
Solucin verdadera

0.22
0.2
0.18
0.16
0.14
0.12
0.4

0.45

0.5

0.55
Tiempo t (seg)

Seccin de Estudios de Posgrado e Investigacin

14

0.6

0.65

0.7

Escuela Superior de Ingeniera Mecnica y Elctrica.

Sistemas Lineales.
Comparacin del Mtodo de Heun con la solucin analtica
h=2-2

0.2213

h=2-4
h=2-6

0.2213

h=2-8

x(t)

0.2213

h=2-10
h=2-12
Solucin verdadera

0.2213
0.2213
0.2213
0.2213
0.2213

0.6049

0.6049

0.6049
0.6049
Tiempo t (seg)

0.6049

0.6049

B) SOLUCIN POR EL MTODO TRAPEZOIDAL

Se usa exactamente el mismo cdigo de MatLab para graficar que en el
mtodo de Heun mostrado anteriormente, donde t1 a t6 y x1 a x6 son las
salidas

del

programa

en

ForTran,

h=22 , 24 , 26 ,28 , 210 y 212

para

respectivamente. Se observa que

Comparacin del Mtodo Trapezoidal con la solucin analtica

1
0.9

h=2-2

0.8

h=2-6

h=2-4
h=2-8

0.7

h=2-10

x(t)

0.6

h=2-12
Solucin verdadera

0.5
0.4
0.3
0.2
0.1
0

0.1

0.2

0.3

0.4
0.5
0.6
Tiempo t (seg)

Seccin de Estudios de Posgrado e Investigacin

15

0.7

0.8

0.9

Escuela Superior de Ingeniera Mecnica y Elctrica.

Sistemas Lineales.
Se puede observar que el mtodo de Heun es ms exacto en ste caso, ya que
es un sistema no rgido y el mtodo de Heun tiene una mejor aproximacin de
la pendiente.
Comparacin del Mtodo Trapezoidal con la solucin analtica
0.395
0.39
0.385

x(t)

0.38
0.375

h=2-2

0.37

h=2-4

0.365

h=2-8

h=2-6
h=2-10

0.36
0.355
0.71

h=2-12
Solucin verdadera
0.72

0.73

0.74
0.75
Tiempo t (seg)

0.76

Seccin de Estudios de Posgrado e Investigacin

16

0.77

0.78

Escuela Superior de Ingeniera Mecnica y Elctrica.

Sistemas Lineales.

EJERCICIO 3
x 1=2 x 2+ 2t
x 2=

x 1 ( 0 ) =2

1
x +2t
2 1

x 2 ( 0 ) =0

d e t=0 a t=2 para h=0.05,0 .005 y 0.0005

A) SOLUCIN POR EL MTODO DE EULER

Cdigo usado para graficar en MatLab, donde t1 a t6 y x1 a x6 son las salidas
del programa en ForTran, para

h=0 .05 , 0 . 005 y 0 . 0005

respectivamente. El

cdigo en ForTran del programa para el mtodo de Euler est al final de este
documento (para ver cmo se introducen 2 variables de estado en el cdigo de
ForTran). Si se requiere aadir otro valor de h se modifica la funcin en MatLab,
aadiendo una x7, x8 y t7, t8 con sus etiquetas correspondientes.
function euler_ejercicio3(t1,x1,t2,x2,t3,x3,t4,x4,t5,x5,t6,x6)
%Grfica de la solucin numrica de x1(t) para h=0.05:
plot(t1,x1,'-y')
hold on
%Grfica de la solucin numrica de x2(t) para h=0.05:
plot(t2,x2,'-m')
hold on
%Grfica de la solucin numrica de x1(t) para h=0.005:
plot(t3,x3,'-b')
hold on
%Grfica de la solucin numrica de x2(t) para h=0.005:
plot(t4,x4,'-r')
hold on
%Grfica de la solucin numrica de x1(t) para h=0.0005:
plot(t5,x5,'-g')
hold on
%Grfica de la solucin numrica de x2(t) para h=0.0005:
plot(t6,x6,'-k')
hold on
title('Comparacin las soluciones usando del Mtodo de Euler para el
ejercicio 3');
xlabel('Tiempo t (seg)');
ylabel('x(t)');
legend('x1(t) para h=0.05','x2(t) para h=0.05','x1(t) para
h=0.005','x2(t) para h=0.005','x1(t) para h=0.0005','x2(t) para
h=0.0005');
end

Seccin de Estudios de Posgrado e Investigacin

17

Escuela Superior de Ingeniera Mecnica y Elctrica.

Sistemas Lineales.
4

Comparacin las soluciones usando del Mtodo de Euler para el ejercicio 3

x1(t) para h=0.05

x2(t) para h=0.05
x1(t) para h=0.005
x2(t) para h=0.005
x1(t) para h=0.0005
x2(t) para h=0.0005

x(t)

-1

-2

0.2

0.4

0.6

0.8
1
1.2
Tiempo t (seg)

1.4

1.6

1.8

Se acerca el zoom para observar la aproximacin de cada solucin. Para x1(t)

se tiene la grfica siguiente, como se puede ver, mientras sea menor el paso
tiende a estabilizarse.
Comparacin las soluciones usando del Mtodo de Euler para el ejercicio 3
0.4
x1(t) para h=0.05
x2(t) para h=0.05
x1(t) para h=0.005
x2(t) para h=0.005
x1(t) para h=0.0005
x2(t) para h=0.0005

0.3

x(t)

0.2
0.1
0
-0.1
-0.2
1.45

1.5

1.55
Tiempo t (seg)

Seccin de Estudios de Posgrado e Investigacin

18

1.6

1.65

Escuela Superior de Ingeniera Mecnica y Elctrica.

Sistemas Lineales.
Para la solucin x2(t) se tiene
Comparacin las soluciones usando del Mtodo de Euler para el ejercicio 3
-0.05
x1(t) para h=0.05
x2(t) para h=0.05
x1(t) para h=0.005
x2(t) para h=0.005
x1(t) para h=0.0005
x2(t) para h=0.0005

-0.1

x(t)

-0.15

-0.2

-0.25

0.45

0.5

0.55

0.6

0.65
0.7
Tiempo t (seg)

0.75

0.8

0.85

B) SOLUCIN POR EL MTODO TRAPEZOIDAL

Se usa un cdigo en MatLab similar al que se us en la solucin del mtodo de
Euler. Y se obtienen las grficas siguientes:
Comparacin las soluciones usando del Mtodo Trapezoidal para el ejercicio 3
4
x1(t) para h=0.05
x2(t) para h=0.05
x1(t) para h=0.005
x2(t) para h=0.005
x1(t) para h=0.0005
x2(t) para h=0.0005

x(t)

-1

-2

0.2

0.4

0.6

0.8
1
1.2
Tiempo t (seg)

Seccin de Estudios de Posgrado e Investigacin

19

1.4

1.6

1.8

Escuela Superior de Ingeniera Mecnica y Elctrica.

Sistemas Lineales.

Haciendo zoom en la solucin x1(t) se observa una mayor exactitud que en el

mtodo de Euler
Comparacin las soluciones usando del Mtodo Trapezoidal para el ejercicio 3
x1(t) para h=0.05
x2(t) para h=0.05
x1(t) para h=0.005
x2(t) para h=0.005
x1(t) para h=0.0005
x2(t) para h=0.0005

0.1242
0.124

x(t)

0.1238
0.1236
0.1234
0.1232
0.123
0.1228
1.632

1.6322

1.6324

1.6326
1.6328
Tiempo t (seg)

1.633

1.6332

Mientras que para la solucin x2(t) se obtiene la misma exactitud que en la

solucin x1(t)
Comparacin las soluciones usando del Mtodo Trapezoidal para el ejercicio 3
1.3278
1.3276
1.3274

x1(t) para h=0.05

x2(t) para h=0.05
x1(t) para h=0.005
x2(t) para h=0.005
x1(t) para h=0.0005
x2(t) para h=0.0005

x(t)

1.3272
1.327
1.3268
1.3266
1.3264
1.5245 1.5246 1.5247 1.5248 1.5249 1.525
Tiempo t (seg)

Seccin de Estudios de Posgrado e Investigacin

20

1.5251 1.5252

Escuela Superior de Ingeniera Mecnica y Elctrica.

Sistemas Lineales.

C) SOLUCIN POR EL MTODO DE RUNGE KUTTA DE 4TO

ORDEN
Usando la misma funcin para graficar en MatLab se obtiene la siguiente
grfica.
Comparacin las soluciones usando del Mtodo de Runge Kutta de 4to Orden para el ejercicio 3
4
x1(t) para h=0.05
x2(t) para h=0.05
x1(t) para h=0.005
x2(t) para h=0.005
x1(t) para h=0.0005
x2(t) para h=0.0005

x(t)

-1

-2

0.2

0.4

0.6

0.8
1
1.2
Tiempo t (seg)

1.4

1.6

1.8

Haciendo zoom en las grficas de x1(t) y x2(t) se observa que el mtodo de

Runge Kutta es un poco menos preciso que el mtodo Trapezoidal, con
resultados muy similares, pero mucho ms preciso que el mtodo de Euler para
este tipo de sistema.

Seccin de Estudios de Posgrado e Investigacin

21

Escuela Superior de Ingeniera Mecnica y Elctrica.

Sistemas Lineales.
Comparacin las soluciones usando del Mtodo de Runge Kutta de 4to Orden para el ejercicio 3
0.162
x1(t) para h=0.05
0.16
x2(t) para h=0.05
x1(t) para h=0.005
0.158
x2(t) para h=0.005
x1(t) para h=0.0005
0.156
x2(t) para h=0.0005

x(t)

0.154
0.152
0.15
0.148
0.146
0.144
0.142
1.636

1.638

1.64
1.642
Tiempo t (seg)

1.644

1.646

1.648

Comparacin las soluciones usando del Mtodo de Runge Kutta de 4to Orden para el ejercicio 3

0.92

x(t)

0.915

x1(t) para h=0.05

x2(t) para h=0.05
x1(t) para h=0.005
x2(t) para h=0.005
x1(t) para h=0.0005
x2(t) para h=0.0005

0.91

0.905

0.9
1.374

1.375

1.376

1.377
1.378
1.379
Tiempo t (seg)

1.38

1.381

D) SOLUCIN POR EL MTODO DE GEAR DE 1ER ORDEN O

EULER REGRESIVO
Usando la misma funcin para graficar en MatLab se obtiene la siguiente
figura.

Seccin de Estudios de Posgrado e Investigacin

22

Escuela Superior de Ingeniera Mecnica y Elctrica.

Sistemas Lineales.
Comparacin las soluciones usando del Mtodo de Gear de 1er Orden para el ejercicio 3
4
x1(t) para h=0.05
x2(t) para h=0.05
x1(t) para h=0.005
x2(t) para h=0.005
x1(t) para h=0.0005
x2(t) para h=0.0005

x(t)

-1

-2

0.2

0.4

0.6

0.8
1
1.2
Tiempo t (seg)

1.4

1.6

1.8

Haciendo zoom en ambas soluciones se tiene que es menos preciso que el

mtodo Trapezoidal y que el Runge Kutta de 4to orden, pero es mejor que el
Euler explcito.
Comparacin las soluciones usando del Mtodo de Gear de 1er Orden para el ejercicio 3
x1(t) para h=0.05
x2(t) para h=0.05
x1(t) para h=0.005
x2(t) para h=0.005
x1(t) para h=0.0005
x2(t) para h=0.0005

0.62

0.615

x(t)

0.61
0.605

0.6

0.595

1.18

1.19

1.2

1.21
1.22
1.23
Tiempo t (seg)

Seccin de Estudios de Posgrado e Investigacin

23

1.24

1.25

1.26

Escuela Superior de Ingeniera Mecnica y Elctrica.

Sistemas Lineales.
Comparacin las soluciones usando del Mtodo de Gear de 1er Orden para el ejercicio 3
0.2
x1(t) para h=0.05
x2(t) para h=0.05
0.15
x1(t) para h=0.005
x2(t) para h=0.005
x1(t) para h=0.0005
0.1
x2(t) para h=0.0005

x(t)

0.05
0
-0.05
-0.1
1.58 1.585 1.59 1.595

1.6 1.605 1.61 1.615 1.62 1.625 1.63

Tiempo t (seg)

EJERCICIO 4
x 1=2 x 1+ x2 +100
x 2=104 x 1104 x 2+ 50

x 1 ( 0 ) =0
x 2 ( 0 ) =0

d e t=0 a t =10 para h=0.0002

A) DETERMINAR h MXIMO CON EL MTODO DE GEAR DE 1ER
ORDEN
La siguiente grfica se obtuvo usando pasos de h=0.0001,0.0002 y 0.01

Seccin de Estudios de Posgrado e Investigacin

24

Escuela Superior de Ingeniera Mecnica y Elctrica.

Sistemas Lineales.
Comparacin las soluciones usando del Mtodo de Gear de 1er Orden para el ejercicio 4
120

100

x(t)

80

x1(t) para h=0.0001

x2(t) para h=0.0001
x1(t) para h=0.0002
x2(t) para h=0.0002
x1(t) para h=0.01
x2(t) para h=0.01

60

40

20

4
5
6
Tiempo t (seg)

10

Haciendo un zoom en la grfica se puede apreciar que tan estable es para los
valores de h dados, se puede deducir que para h=0.01 tiene ms error en la
aproximacin. Se selecciona el paso mximo como h=0.0002
Comparacin las soluciones usando del Mtodo de Gear de 1er Orden para el ejercicio 4

93.5
93.48

x(t)

93.46
x1(t) para h=0.0001
x2(t) para h=0.0001
x1(t) para h=0.0002
x2(t) para h=0.0002
x1(t) para h=0.01
x2(t) para h=0.01

93.44
93.42
93.4
93.38
2.7285

2.7285

2.7286

2.7287
2.7287
Tiempo t (seg)

2.7287

Seccin de Estudios de Posgrado e Investigacin

25

2.7288

2.7289

Escuela Superior de Ingeniera Mecnica y Elctrica.

Sistemas Lineales.
Comparacin las soluciones usando del Mtodo de Gear de 1er Orden para el ejercicio 4
42.195
42.19
42.185
x1(t) para h=0.0001
x2(t) para h=0.0001
x1(t) para h=0.0002
x2(t) para h=0.0002
x1(t) para h=0.01
x2(t) para h=0.01

x(t)

42.18
42.175
42.17
42.165
42.16
0.5475

0.548

0.5485

0.549 0.5495
Tiempo t (seg)

0.55

0.5505

0.551

B Y C) SOLUCIN CON EL MTODO DE EULER PARA 2hmax y

hmax/2
Claramente se aprecia en la grfica obtenida que para un valor mayor de la h
mxima determinada por el mtodo de Gear en el inciso anterior el sistema de
vuelve inestable, mientras que para valores menores a h mxima el sistema no
presenta inestabilidad.
2

6
Comparacin
las soluciones usando del Mtodo de Euler para el ejercicio 4
x 10

x1(t) para h=0.0001

x2(t) para h=0.0001
x1(t) para h=0.0004
x2(t) para h=0.0004

1.5

x(t)

0.5

-0.5

-1

4
5
6
Tiempo t (seg)

10

Haciendo un zoom en ambas grficas, se puede ver ms claramente esto. Para

h=0.0004 se comporta de la forma siguiente

Seccin de Estudios de Posgrado e Investigacin

26

Escuela Superior de Ingeniera Mecnica y Elctrica.

Sistemas Lineales.
Comparacin las soluciones usando del Mtodo de Euler para el ejercicio 4
6000

x1(t) para h=0.0001

x2(t) para h=0.0001
x1(t) para h=0.0004
x2(t) para h=0.0004

4000

x(t)

2000

-2000

-4000

5
Tiempo t (seg)

8
-3

x 10

Y para h=0.0001 de sta otra forma:

Comparacin las soluciones usando del Mtodo de Euler para el ejercicio 4
x1(t) para h=0.0001
x2(t) para h=0.0001
x1(t) para h=0.0004
x2(t) para h=0.0004

0.63
0.62
0.61

x(t)

0.6
0.59
0.58
0.57
0.56
0.55
6.085

6.09

6.095
6.1
Tiempo t (seg)

6.105

Seccin de Estudios de Posgrado e Investigacin

27

6.11

6.115
-3

x 10

Escuela Superior de Ingeniera Mecnica y Elctrica.

Sistemas Lineales.

CDIGO DE LOS MTODOS EN FORTRAN 90

!****************************************************************************
!
! PROGRAMA: euler (Adams - Bashford 1er orden)
!
! OBJETIVO: Hallar el valor numrico aproximado para la solucin de una
! ecuacin diferencial de la forma y'=f(x,y)
!
!****************************************************************************
program datos
implicit none
double precision, dimension(100000) :: x,t
real h,x0,df,t0,t1
integer i,puntos
print *, '****************************************'
print *, '
MTODO
DE EULER
'
print *, '****************************************'
print *, 'Tiempo inicial: '
read *, t0
print *, 'Tiempo final:'
read *, t1
print *, 'Valor inicial de x(t): '
read *, x0
print *, 'Valor de h: '
read *, h
10 format (T8, 't', T22 , 'x(t)' )
20 format (T2, F10.6, T16, F10.6)
puntos = (t1 - t0)/h
x(1) = x0
t(1) = t0
do i=1,puntos
call Yprima(x,i,df)
call Euler(x,i,h,t,df)
end do
open(7,FILE='EULER_log.TXT')
write(7,10)
do i=1,puntos+1
write(7,20) t(i),x(i)
end do
close(7)
end program datos
subroutine Yprima(x,i,df)
double precision, dimension(100000) :: x
df = x(i)**2
return
end subroutine

Seccin de Estudios de Posgrado e Investigacin

28

Escuela Superior de Ingeniera Mecnica y Elctrica.

Sistemas Lineales.
subroutine Euler(x,i,h,t,df)
double precision, dimension(100000) :: x,t
x(i+1) = x(i) + df*h
t(i+1) = t(i) + h
return
end subroutine
!****************************************************************************
!
! PROGRAMA: Trapezoidal (Adams - Moulton Implcito)
!
! OBJETIVO: Hallar el valor numrico aproximado para la solucin de una
! ecuacin diferencial de la forma y'=f(x,y)
!
!****************************************************************************
program trapezoidal
implicit none
double precision, dimension(100000) :: x,t
double precision h,x0,df,t0,t1
integer i,ciclos
print *, '****************************************'
print *, '
METODO
TRAPEZOIDAL
'
print *, '****************************************'
print *, 'tiempo inicial: '
read *, t0
print *, 'tiempo final:'
read *, t1
print *, 'valor inicial de x(t): '
read *, x0
print *, 'valor de h: '
read *, h
10 format (t8, 't', t22 , 'x(t)' )
20 format (t2, f10.6, t16, f12.6)
ciclos = (t1 - t0)/h
x(1) = x0
t(1) = t0
do i=1,ciclos
call procesotrapecio(x,t,i,h)
end do
open(7,file='trapecio_log.txt')
write(7,10)
do i=1,ciclos+1
write(7,20) t(i),x(i)
end do
close(7)
end program trapezoidal

Seccin de Estudios de Posgrado e Investigacin

29

Escuela Superior de Ingeniera Mecnica y Elctrica.

Sistemas Lineales.
subroutine procesotrapecio(x,t,i,h)
double precision, dimension(100000) :: x,t
double precision alpha,xn,h
alpha=(x(i)**2)+(2/h)*x(i)
call newton(x,i,h,xn,alpha)
x(i+1) = xn
t(i+1) = t(i) + h
return
end subroutine
subroutine newton(x,i,h,xn,alpha)
double precision, dimension(100000) :: x
double precision xn,fx,dfx,alpha,h
real :: tol = 0.001
integer :: nit = 100
integer it
it = 1
xn = x(i)
do
fx = xn**2 - (2/h)*xn + alpha
dfx = 2*xn - (2/h)
xn = xn - (1/dfx)*(fx)
fx = xn**2 - (2/h)*xn + ((x(i)**2) + x(i)*(2/h))
if(abs(fx)<=tol) exit
if(it>=nit)exit
it = it + 1
end do
return
end
!****************************************************************************
!
! PROGRAMA: runge kutta de 4to orden
!
! OBJETIVO: hallar el valor numrico aproximado para la solucin de una
! ecuacin diferencial de la forma y'=f(x,y)
!
!****************************************************************************
program rungekutta
implicit none
double precision, dimension(100000) :: x,t
real h,x0,t0,tf
double precision k1,k2,k3,k4
integer i,nc
print *, '****************************************'
print *, ' METODO
RUNGE KUTTA 4TO ORDEN
print *, '****************************************'

'

print *, 'introduce el tiempo inicial: '

Seccin de Estudios de Posgrado e Investigacin

30

Escuela Superior de Ingeniera Mecnica y Elctrica.

Sistemas Lineales.
read *, t0
print *, 'introduce el tiempo final: '
read *, tf
print *, 'introduce el valor inicial de x: '
read *, x0
print *, 'introduce el valor del paso: '
read *, h
10 format (t8, 't', t22 , 'x'/ t6, 5('-'), t20, 5('-') /)
20 format (t2, f12.6, t16, f12.6)
nc = (tf - t0)/h
x(1) = x0
t(1) = t0
do i=1,nc
call yprima(x,i,h,k1,k2,k3,k4)
call procesorungek(x,i,h,t,k1,k2,k3,k4)
end do
open(7,file='kutta4_log.txt')
write(7,10)
do i=1,nc+1
write(7,20) t(i),x(i)
end do
close(7)
end program
subroutine yprima(x,i,h,k1,k2,k3,k4)
double precision, dimension(100000) :: x
double precision k1,k2,k3,k4
k1 = x(i)**2
k2 = (x(i) + (h/2)*k1)**2
k3 = (x(i) + (h/2)*k2)**2
k4 = (x(i) + h*k3)**2
return
end
subroutine procesorungek(x,i,h,t,k1,k2,k3,k4)
double precision, dimension(100000) :: x,t
double precision k1,k2,k3,k4
x(i+1) = x(i) + (h/6)*(k1 + 2*k2 + 2*k3 + k4)
t(i+1) = t(i) + h
return
end

!****************************************************************************
!
! PROGRAMA: Gear de 1er orden o Euler regresivo
!
! OBJETIVO: hallar el valor numrico aproximado para la solucin de una
! ecuacin diferencial de la forma y'=f(x,y)
!
!****************************************************************************
program trapezoidal
implicit none

Seccin de Estudios de Posgrado e Investigacin

31

Escuela Superior de Ingeniera Mecnica y Elctrica.

Sistemas Lineales.
double precision, dimension(100000) :: x,t
double precision h,x0,df,t0,t1
integer i,ciclos
print *, '****************************************'
print *, '
MTODO
DE
GEAR
'
print *, '****************************************'
print *, 'tiempo inicial: '
read *, t0
print *, 'tiempo final:'
read *, t1
print *, 'valor inicial de x(t): '
read *, x0
print *, 'valor de h: '
read *, h
10 format (t8, 't', t22 , 'x(t)' )
20 format (t2, f10.6, t16, f12.6)
ciclos = (t1 - t0)/h
x(1) = x0
t(1) = t0
do i=1,ciclos
call procesotrapecio(x,t,i,h)
end do
open(7,file='gear_log.txt')
write(7,10)
do i=1,ciclos+1
write(7,20) t(i),x(i)
end do
close(7)
end program trapezoidal
subroutine procesotrapecio(x,t,i,h)
double precision, dimension(100000) :: x,t
double precision alpha,xn,h
alpha=(1/h)*x(i)
call newton(x,i,h,xn,alpha)
x(i+1) = xn
t(i+1) = t(i) + h
return
end subroutine
subroutine newton(x,i,h,xn,alpha)
double precision, dimension(100000) :: x
double precision xn,fx,dfx,alpha,h
real :: tol = 0.001
integer :: nit = 100
integer it

Seccin de Estudios de Posgrado e Investigacin

32

Escuela Superior de Ingeniera Mecnica y Elctrica.

Sistemas Lineales.
it = 1
xn = x(i)
do
fx = xn**2 - (1/h)*xn + alpha
dfx = 2*xn - (1/h)
xn = xn - (1/dfx)*(fx)
fx = xn**2 - (1/h)*xn + alpha
if(abs(fx)<=tol) exit
if(it>=nit)exit
it = it + 1
end do
return
end

!****************************************************************************
!
! PROGRAMA: Heun
!
! OBJETIVO: hallar el valor numrico aproximado para la solucin de una
! ecuacin diferencial de la forma y'=f(x,y)
!
!****************************************************************************

program heun
implicit none
double precision, dimension(100000) :: x,t
double precision yprima,k1,k2,aux
real h,x0,df,t0,t1
integer i,ciclos

print *, '****************************************'
print *, '
metodo
de heun
'
print *, '****************************************'
print *, 'tiempo inicial: '
read *, t0
print *, 'tiempo final:'
read *, t1
print *, 'valor inicial de x(t): '
read *, x0
print *, 'valor de h: '
read *, h
10 format (t8, 't', t22 , 'x(t)' )
20 format (t2, f10.6, t16, f10.6)
ciclos = (t1 - t0)/h
x(1) = x0

Seccin de Estudios de Posgrado e Investigacin

33

Escuela Superior de Ingeniera Mecnica y Elctrica.

Sistemas Lineales.
t(1) = t0
do i=1,ciclos
t(i+1) = t(i) + h
k1=h*yprima(t(i),x(i))
aux=x(i)+k1
k2=h*yprima(t(i+1),aux)
x(i+1) = x(i)+.5*(k1+k2)
end do
open(7,file='heun_log.txt')
write(7,10)
do i=1,ciclos+1
write(7,20) t(i),x(i)
end do
close(7)
end program heun
double precision function yprima(tt,xx)
double precision tt,xx
yprima=xx-tt**3+3*tt**2
return
end function yprima

!****************************************************************************
!
! PROGRAM: Euler
!
! PURPOSE: Ejemplo para introducir datos cuando se tiene un sistema de
2 variables de estado.
!
!****************************************************************************
program euler
implicit none
double precision, dimension(100000) :: x1,x2,t
real h,x01,x02,df1,df2,t0,t1
integer i,ciclos
print *, '****************************************'
print *, '
METODO
DE EULER
'
print *, '****************************************'
print *, 'tiempo inicial: '
read *, t0
print *, 'tiempo final:'
read *, t1
print *, 'valor inicial de x1(t): '
read *, x01
print *, 'valor inicial de x2(t): '
read *, x02
print *, 'valor de h: '
read *, h

Seccin de Estudios de Posgrado e Investigacin

34

Escuela Superior de Ingeniera Mecnica y Elctrica.

Sistemas Lineales.
10 format (t8, 't', t22 , 'x1(t)', t36, 'x2(t)')
20 format (t2, f14.6, t16, f14.6, t30, f14.6)
ciclos = (t1 - t0)/h
x1(1) = x01
x2(1) = x02
t(1) = t0
do i=1,ciclos
call yprima(x1,x2,t,i,df1,df2)
call procesoeuler(x1,x2,i,h,t,df1,df2)
end do
open(7,file='euler_log.txt')
write(7,10)
do i=1,ciclos+1
write(7,20) t(i),x1(i),x2(i)
end do
close(7)
end program euler
subroutine yprima(x1,x2,t,i,df1,df2)
double precision, dimension(100000) :: x1,x2,t
df1 = -2*x2(i)+2*t(i)**2
df2 = .5*x1(i)+2*t(i)
return
end subroutine
subroutine procesoeuler(x1,x2,i,h,t,df1,df2)
double precision, dimension(100000) :: x1,x2,t
t(i+1) = t(i) + h
x1(i+1) = x1(i) + df1*h
x2(i+1) = x2(i) + df2*h
return
end subroutine

CONLCUSIONES
Como pudo observarse en el transcurso del desarrollo de ste documento es
que no siempre una aproximacin por un mtodo implcito para un sistema no
rgido es la ms exacta, de ello depende de qu forma ste fragmentada la
pendiente.
Se puede concluir tambin que en los sistemas rgidos como es el caso del
ejercicio 4, se deben hacer las pruebas suficientes para determinar su
estabilidad, para encontrar el h mximo o los intervalos de tiempo, no siempre
es recomendable hacerlo a prueba y error cmo se hizo en ste documento, se
debe tener un conocimiento ms a fondo de esto, y se llega a que la rigidez de
un sistema depende de los valores propios (modos), lo que probablemente se
examine a ms detalle en trabajos posteriores.

Seccin de Estudios de Posgrado e Investigacin

35

Escuela Superior de Ingeniera Mecnica y Elctrica.

Sistemas Lineales.

BIBLIOGRAFA
[1] Computer-Aided Analysis of Electronic Circuits, Algorithms and
Computational Techniques, by Leon O. Ch. (1975).
[2] Numerical Methods for Engineers by Steven Chapra and Raymond
Canale (Jan 24, 2014).
[3] Introduction to Programming with Fortran: With Coverage of Fortran 90, 95,
2003, 2008 and 77 by ian chivers and Jane Sleightholme (Feb 9, 2012).

Seccin de Estudios de Posgrado e Investigacin

36