Araiza_Rincón_Sánchez_T8 - Axel

Vista previa del material en texto

Análisis Numérico
Grupo: 13
UNAM / FI / DCB(DIE)
Semestre 2021-1
Tarea 8
Ecuaciones Diferenciales y Sistemas de Ecuaciones Diferenciales
Araiza Alfredo, Rincón Carlos, Sánchez Axel
1. Gráfica comparativa de la solución de la ecuación diferen-
cial de primer orden: xy′ + y = ey/x
Euler y Euler mejorado con h = 0.1
Runge-Kutta 4to orden con h = 0.3
1 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9
1
1.1
1.2
1.3
1.4
1.5
1.6
1.7
1.8
1.9
2
2.1
2.2
Euler
Euler mejorado
RK4 clásico
AN21-1x13: Tarea 8 Araiza Alfredo, Rincón Carlos, Sánchez Axel
2. Gráfica comparativa de la solución de la ecuación diferen-
cial de segundo orden: y′′ + 3y′ + 2y = x + 1
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
0
0.025
0.05
0.075
0.1
0.125
0.15
0.175
0.2
0.225
0.25
0.275
0.3
Euler
Heun
RK4
2
AN21-1x13: Tarea 8 Araiza Alfredo, Rincón Carlos, Sánchez Axel
3. Gráfica comparativa de la solución de la ecuación diferen-
cial con valores en la frontera: 2u′′ + u′ + 3u = 23sin(
t
2)
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
X
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
1.1
1.2
1.3
1.4
Y
Límites
Primera adivinanza
Segunda adivinanza
Metodo de Disparo
Diferencias Finitas
4. Descripción de trabajo en equipo
Para obtener los resultados, nos repartimos equitativamente todo el trabajo, Alfredo se encargó
de realizar y modificar los programas para el inciso 2, es decir, todos aquellos que se ocupan para
poder realizar los distintos métodos para un sistema de ecuaciones diferenciales, Carlos realizó los
programas para poder efectuar correctamente los métodos que nos permiten obtener la solución
de una ecuación diferencial de primer orden, además del desarrollo a mano y pasar la información
a LATEXy Axel se encargó de realizar todos los programas con respecto al método de disparo y
diferencias finitas. Al final realizamos una videollamada para verificar que el archivo y su contenido
sean los correctos y para aclarar las dudas que surgieron en el proceso.
3
AN21-1x13: Tarea 8 Araiza Alfredo, Rincón Carlos, Sánchez Axel
5. Desarrollo
4
AN21-1x13: Tarea 8 Araiza Alfredo, Rincón Carlos, Sánchez Axel
5
AN21-1x13: Tarea 8 Araiza Alfredo, Rincón Carlos, Sánchez Axel
6. Códigos utilizados
’Código Principal 1’
1 clc % Borra Command Window
2 clear all % Borra Workspace
3 close all % Cierra ventanas
4
5 f1 = @(x,y) (exp(y/x)-y)/(x); % Ecuación a resolver
6 h1 = 0.1; % Paso para Euler y Euler mejorado
7 h2 = 0.3; % Paso para RK4
8 x = 1:h1:1.9; % Dominio
9 x2 = 1:h2:1.9; % Dominio
10 y0 = 1; % Condición inicial
11
12 %% Solución por el método de Euler
13 y = eulerEDO(f1,x,y0)
14 plot(x,y,’Color’,’c’,’MarkerSize’,6,’LineWidth’,2)
15 %% Solución por el método de Euler mejorado
16 y1 = eulermejoradoEDO(f1,x,y0)
17 hold on
18 plot(x,y1,’Color’,’g’,’MarkerSize’,6,’LineWidth’,2)
19 %% Solución por el método de RK4 clásico
20 y2 = RK4clasico(f1,x2,y0)
21 hold on
22 plot(x2,y2,’s’,’Color’,’m’,’MarkerSize’,6,’LineWidth’,2)
23 xlabel(’x’,’interpreter’,’latex’,’Fontsize’,16)
24 ylabel(’y’,’interpreter’,’latex’,’Fontsize’,16)
25 legend(’Euler’,’Euler mejorado’,’RK4 clásico’)
’Código Euler’
1 function y = eulerEDO(f,x,y0)
2 y = 0*x; % Define vector y solución
3 y(1) = y0; % Condición inicial
4 h = x(2) - x(1); % Determinar el paso
5 n = length(x); % Determinar el número de elementos de x
6 for i=1:n-1
7 y(i+1) = y(i) + h*f(x(i),y(i)); % Método de Euler
8 end
9 end
10
’Código Euler mejorado’
1 function y = eulermejoradoEDO(f,x,y0)
2 y = 0*x; % Define vector y solución
6
AN21-1x13: Tarea 8 Araiza Alfredo, Rincón Carlos, Sánchez Axel
3 y(1) = y0; % Condición inicial
4 h = x(2) - x(1); % Determinar el paso
5 n = length(x); % Determinar el número de elementos de x
6 for i=1:n-1
7 k1 = f(x(i),y(i));
8 k2 = f(x(i)+(1/2)*h,y(i)+1/2*k1*h);
9
10 y(i+1) = y(i) + h*k2; % Método de Euler mejorado
11 end
12 end
’Código RK4 clásico’
1 function y = RK4clasico(f,x,y0)
2 y = 0*x; % Vector "y" lleno de ceros (dimensiones de x)
3 y(1) = y0; % Condición inicial
4 h = x(2)-x(1); % Paso
5 n = length(x); % Tamaño del vector "x"
6 for i = 1:n-1
7 k1 = f(x(i),y(i)); % Calcula k1
8 k2 = f(x(i)+(1/2)*h,y(i)+(1/2)*k1*h); % Calcula k2
9 k3 = f(x(i)+(1/2)*h,y(i)+(1/2)*k2*h); % Calcula k3
10 k4 = f(x(i)+h,y(i)+k3*h); % Calcula k4
11
12 y(i+1) = y(i) + (1/6)*h*(k1+2*k2+2*k3+k4); % Calcula y(i+1)
13 end
14 end
’Código Principal 2’
1 clc % Borra Command Window
2 clear all % Borra Workspace
3 close all % Cierra ventanas
4
5 f2 = @(x,u) [ u(2) ; -3*u(2)-2*u(1)+x+1 ]; % Vector con ec. de estado
6 u0 = [0;0]; % Vector con condiciones iniciales
7 x = 0:0.1:1; % Dominio dividido en "n" pasos
8
9 %% Por el método de Euler
10 Euler = EulerSEDO(f2,x,u0)
11 plot(x,Euler(1,:),’Color’,’g’,’MarkerSize’,6,’LineWidth’,2)
12 %% Poel método de Heun
13 Heun = HeunSEDO(f2,x,u0)
14 hold on
15 plot(x,Heun(1,:),’s’,’Color’,’m’,’MarkerSize’,6,’LineWidth’,2)
16 %% Por el método de Runge-Kutta 4 orden
7
AN21-1x13: Tarea 8 Araiza Alfredo, Rincón Carlos, Sánchez Axel
17 RK4 = RK4SEDO(f2,x,u0)
18 hold on
19 plot(x,RK4(1,:),’Color’,’c’,’MarkerSize’,6,’LineWidth’,2)
20 xlabel(’x’,’interpreter’,’latex’,’Fontsize’,16)
21 ylabel(’y’,’interpreter’,’latex’,’Fontsize’,16)
22 legend(’Euler’,’Heun’,’RK4’)
’Código Euler SEDO’
1 function u = EulerSEDO(f,x,u0)
2 u(:,1) = u0; % Primera columna es vector con condiciones iniciales
3 h = x(2) - x(1); % Paso
4 n = length(x); % Tamaño del vector x
5
6 for i = 1:n-1
7 u(:,i+1) = u(:,i) + h*f(x(i),u(:,i)); % Método de Euler SEDO
8 end
9 end
10
’Código Heun SEDO’
1 function u = HeunSEDO(f,x,u0)
2 u(:,1) = u0; % Primera columna es vector con condiciones iniciales
3 h = x(2) - x(1); % Paso
4 n = length(x); % Tamaño del vector x
5
6 for i = 1:n-1
7
8 k1 = f(x(i),u(:,i)); % Calcula k1
9 k2 = f(x(i)+h,u(:,i)+h*k1); % Calcula k2
10
11 u(:,i+1) = u(:,i) + (1/2)*h*(k1+k2); % Método de Heun
12 end
13 end
’Código RK4 SEDO’
1 function u = RK4SEDO(f,x,u0)
2 u(:,1) = u0; % Primera columna es vector con condiciones iniciales
3 h = x(2) - x(1); % Paso
4 n = length(x); % Tamaño del vector x
5
6 for i = 1:n-1
7
8
AN21-1x13: Tarea 8 Araiza Alfredo, Rincón Carlos, Sánchez Axel
8 k1 = f(x(i),u(:,i)); % Calcula k1
9 k2 = f(x(i)+(1/2)*h,u(:,i)+(1/2)*h*k1); % Calcula k2
10 k3 = f(x(i)+(1/2)*h,u(:,i)+(1/2)*h*k2); % Calcula k3
11 k4 = f(x(i)+h,u(:,i)+h*k3); % Calcula k4
12
13 u(:,i+1) = u(:,i) + (1/6)*h*(k1+2*k2+2*k3+k4); % Método RK 4
14 end
15 end
’Código Principal 3’
1 clc % Borra Command Window
2 clear all % Borra Workspace
3 close all % Cierra ventanas
4
5 % Vector con variables de estado
6 f = @(t,x) [ x(2); (1/2)*(-x(2)-3*x(1)+(2/3)*sin(t/2)) ];
7 t = 0:0.05:1; % Dominio de variable independiente
8 % t- variable independiente
9 % x- vector con x1, x2 (variables de estado)
10
11 %% Condiciones de frontera
12 u_0 = 0; % u(0)=0 -> CONOCIDA
13 u_1 = 1;% u(1)=1 -> OBJETIVO
14 figure
15 plot([0 1],[u_0 u_1],’rd’,’MarkerSize’,10,’LineWidth’,2)
16 hold on
17
18 %% Adivinanza 1
19 x0_adv1 = [ u_0 ; 2 ]; % Adivinanza de condiciones iniciales
20
21 x_adv1 = RK4SEDO(f,t,x0_adv1) % Resuelve sistema de EDOs
22 u_1adv1 = x_adv1(1,end);
23 plot(t,x_adv1(1,:),’LineWidth’,2)
24
25 %% Adivinanza 2
26 x0_adv2 = [ u_0 ; 1 ]; % Adivinanza de condiciones iniciales
27
28 x_adv2 = RK4SEDO(f,t,x0_adv2) % Resuelve sistema de EDOs
29 u_1adv2 = x_adv2(1,end);
30 plot(t,x_adv2(1,:),’LineWidth’,2)
31
32 %% Interpolación para x2(0)
33
34 xi = 1; % Valor dónde se interpolará => u(10)=100
35 yi = InterpLagrange([u_1adv1; u_1adv2],[2; 1],xi) % Valor de x2(0)
9
AN21-1x13: Tarea 8 Araiza Alfredo, Rincón Carlos, Sánchez Axel
36
37 %% Resultado Final
38
39 x0_final = [ u_0 ; yi ];
40 x_final = RK4SEDO(f,t,x0_final)
41 plot(t,x_final(1,:),’LineWidth’,2)
42 legend(’Límites’,’Primera adivinanza’,’Segunda adivinanza’,’Método de Disparo’)
43
44 %% Método de diferencias finitas
45
46 A=[-61/32 17/16 0 ; 15/16 -61/32 17/16; 0 15/16 -61/32]; % Vector A
47 B=[0.00259;0.00515;-1.05487]; % Vector B
48
49 u2_3=metCrout(A,B); % Resuelvesistema de ecuaciones lineales
50
51 xDF = 0:0.25:1; % Vector de dominio para método DF
52 solDF = [u_0 u2_3’ u_1] % Vector solución de ED método DF
53
54 plot(xDF,solDF,’pk’,’MarkerSize’,7,’Linewidth’,2)
55 legend(’Límites’,’Primera adivinanza’,’Segunda adivinanza’,’Método de Disparo’,’Diferencias Finitas’)
’Código Interpolación de Lagrange’
1 function yi = InterpLagrange(x,y,xi)
2 % x, y son vectores columna
3 % xi, yi son escalares
4 n1 = length(x); % Determinamos el tamaño del vector (n+1)
5 L = ones(1,n1); % Inicializar la multiplicación
6
7 for i=1:n1 % Para L1,L2,L3...
8 for j=1:n1 % Para el producto
9 if j~= i % Pregunta si j es distinto a i
10 L(i)= L(i)*( xi - x(j) )/( x(i) - x(j) );
11 end
12 end
13 end
14 yi = L*y; % Reemplaza a la sumatoria
15 end
’Código Factorización Crout’
1 function [L,U] = factCrout(A)
2 n= size(A,1); % Tamano de la matriz A
3 L= zeros(n); % Matrizde ceros
4 U= eye(n); % Matriz con 1s (diagonal principal)
5 for i=1: n % i sirve para los renglones en U, columnas en L
10
AN21-1x13: Tarea 8 Araiza Alfredo, Rincón Carlos, Sánchez Axel
6 if i==1 % ¿Primera vuelta?
7 L(1,1) = A(1,1);
8 for j=2:n % j para las columnas de U, renglones de L
9 L(j,1) = A(j,1);
10 U(1,j) = A(1,j)/L(1,1);
11 end
12 else
13 L(i,i) = A(i,i) - L(i,1:i-1)*U(1:i-1,i); %elementos de la diagonal principal
14 for j = i+1 : n
15 L(j,i) = A(j,i) - L(j,1:i-1)*U(1:i-1,i);
16 U(i,j) = (A(i,j) - L(i,1:i-1)*U(1:i-1,j))/L(i,i);
17 end
18 end
19 end
20 L; %matriz L
21 U;%matriz U
22 end
’Código Método Crout’
1 function x = metCrout(A,b)
2 [L,U] = factCrout(A); %Factorizacion [L][U]
3 n = size(A,1);
4
5 % Ly=b
6 y= zeros(n,1); %Vector y
7 y(1)= b(1)/L(1,1);
8 for i= 2 : n
9 y(i) = (b(i) - L(i,1:i-1)*y(1:i-1))/L(i,i);
10 end
11 % Ux=y
12 x= zeros(n,1); % Vector x
13 x(n) = y(n);
14 for i= n-1 : -1 : 1
15 x(i) = y(i)-U(i,i+1:n)*x(i+1:n);
16 end
17 end
*Se usó el mismo programa de RK4SEDO para los principales 2 y 3.
11