Escuela Superior Politécnica del Litoral

Facultad de Ciencias Naturales y Matemáticas

Departamento de Matemáticas

Análisis NuméricoProyecto de curso

Nombre del profesor: M.Sc. Luis Rodríguez Ojeda

Nombre de los integrantes:

Diego Abel Jacho Alvarado

Javier Steveen Carriel Pin

Paralelo: 1

Año lectivo: 2015 - 2016

Objetivos del proyecto

a) Aplicar los conocimientos y técnicas vistas en clase y en el laboratorio para obtener una solución computacional aceptable para problemas propuestos

b) Desarrollar experiencia en aplicaciones computacionales interactuando en equipos

Tema 1Calcular aproximadamente el área de corte de la manzana del siguiente gráfico:

Primero le agregamos ejes coordenados para que sea más visible la obtención de los puntos de interés

X

y

Del gráfico, de manera aproximada tomamos los siguientes 16 puntos:

(11.5, 16), (13,17.1), (15.1,17.9), (20,16), (22.6,12),(22.8,7),(20,3),(18,1),(15.6,1.6),(14,0.8),(11,0.3),(4.5,4)

(2.8, 11), (6,16), (10,17),(13.5,16)

Y se le asignan 16 valores para el parámetro t

t: 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15,16

Para graficar los trazadores cúbicos y graficar el área respectiva se toma 100 puntos del parámetro t, subdividiéndolo

Para realizar la práctica computacional se desarrolló el siguiente algoritmo en lenguaje Matlab

function [P]=trazador(x,y,k)% Trazador Cúbico Natural% p = d(i)(x-x(i))^3+c(i)(x-x(i))^2+b(i)(x-x(i))+a(i), n>2% Entrega k valores de la función a que se intenta aproximar% k puntos del trazadorn=length(x);if n<3 p=[ ]; returnendfor i=1:n-1 h(i)=x(i+1)-x(i);endfor i=1:n-1 m(i)=(y(i+1)-y(i))./h(i);endfor i=1:n-1 a(i)=y(i);enddo=zeros(n,1);do(1)=1; do(n)=1;for i=2:n-1 do(i)=2.*(h(i-1)+h(i));enddn=zeros(n,1);for i=3:n dn(i)=h(i-1);enddm=zeros(n,1);for i=1:n-2 dm(i)=h(i);endB=[dm do dn];b=[-1 0 1];MC=spdiags(B,b,n,n);R=zeros(n,1);R(1)=0; R(n)=0;for i=2:n-1 R(i)=3.*(m(i)-m(i-1));endc=inv(MC)*R;for i=1:n-1

b(i)=m(i)-(h(i)./3).*(2.*c(i)+c(i+1));endfor i=1:n-1 d(i)=(c(i+1)-c(i))./(3.*h(i));ende=(x(n)-x(1))./(k-1);g1=[x(1):e:x(n)];u=length(g1);j=1;for i=1:u-1 if g1(i)<x(j+1) P(i)=a(j)+b(j).*(g1(i)-x(j))+c(j).*(g1(i)-x(j)).^2+d(j).*(g1(i)-x(j)).^3; else j=j+1; P(i)=a(j)+b(j).*(g1(i)-x(j))+c(j).*(g1(i)-x(j)).^2+d(j).*(g1(i)-x(j)).^3; endendi=u; P(i)=a(j)+b(j).*(g1(i)-x(j))+c(j).*(g1(i)-x(j)).^2+d(j).*(g1(i)-x(j)).^3;end

Y en la ventana de comandos de Matlab se realizó el respectivo procedimiento usando los datos anteriormente mencionados para hacer y graficar los trazadores y a su vez graficar los puntos escogidos.

>> t=[0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15];

>> x=[11.5 13 15.1 20 22.6 22.8 20 16 13.6 12 09 2.5 0.8 4 8 11.5];

>> y=[16 17.1 17.9 16 12 7 3 1 1.6 0.8 0.3 4 11 16 17 16];

>> [X]=trazador(t,x,100);

>> [Y]=trazador(t,y,100);

>> figure();

>>plot(x,y,'o','Color',[0,0,0],'MarkerEdgeColor','k','MarkerFaceColor','r','MarkerSize',5),hold on, grid on;

>> hold on, grid on, plot(X,Y);

>>

De lo cual se obtuvo el siguiente gráfico

Hemos logrado obtener un gráfico aproximado a la manzana. Ahora necesitamos calcular su área, para aquello utilizamos la función dada en clase agauss

function [S]=agauss(x,y)n=length(x);S=0;for i=1:n-1 S=S+x(i).*y(i+1);endS=S+x(n).*y(1);for i=1:n-1 S=S-x(i+1).*y(i);endS=S-x(1).*y(n);S=abs(0.5.*S);

Tomando los trazadores X,Y que ya utilizamos para graficar, los usamos para hallar el área deseada:

>> [S]=agauss(X,Y)

S =

311.2605

Obtenemos que el área es aproximadamente 311.2605 unidades cuadradas.

Primero utilizando los datos que el problema proporciona, tenemos:

x ' '+4 (x2−1 ) x'+x=0

Luego se realiza el cambio de variable respectivo para resolver la EDO:

z=x '

Entonces:

z '=x ' '

Además las condiciones iniciales pasarían a ser:

x (0 )=0.75 z (0 )=0

Finalmente tenemos el sistema de EDO’s

x '=z

z '+4 (x2−1 ) z+x=0

Despejando tenemos:

x '=z

z '=4 (1−x2 ) z−x

Para resolverlo usamos el siguiente algoritmo desarrollado en Matlab implementando el método Runge-Kutta de cuarto orden:

function [u,v,w]=rk42(f, g, x, y, z, h, m)for i=1:m k1y=h*f(x,y,z); k1z=h*g(x,y,z); k2y=h*f(x+h/2,y+k1y/2,z+k1z/2); k2z=h*g(x+h/2,y+k1y/2,z+k1z/2); k3y=h*f(x+h/2,y+k2y/2,z+k2z/2); k3z=h*g(x+h/2,y+k2y/2,z+k2z/2); k4y=h*f(x+h,y+k3y,z+k3z); k4z=h*g(x+h,y+k3y,z+k3z); y=y+1/6*(k1y+2*k2y+2*k3y+k4y);

z=z+1/6*(k1z+2*k2z+2*k3z+k4z); x=x+h; u(i)=x; v(i)=y; w(i)=z;endfigure()plot(u,v,'o','Color',[0,0,0],'MarkerEdgeColor','k','MarkerFaceColor','r','MarkerSize',5),hold on, grid on;axis([ 0 1 0 1])end

Este algoritmo requiere que se ingrese funciones f y g, las cuales representan a las derivadas del sistema de forma explícita, también se requieren los puntos iniciales; t para la variable independiente, x para la variable x(t) y z para la variable z(t).

Además el algoritmo requiere ingresar un valor de h que corresponde a la distancia entre los puntos de la variable independiente t y m es el número de puntos correspondientes a la solución deseada.

Las salidas de este algoritmo son los vectores, “u”, ”v” y “w”; el vector u contiene todos los valores de la variable independiente t, el vector v contiene los valores de x(t) y el vector w contiene los valores de z(t).

El algoritmo se encarga de graficar el conjunto de puntos de la solución que son formados por los vectores “v” y “t”

Para utilizar el algoritmo mencionado, definimos las funciones f y g:

f ( t , x , z )=zg ( t , x , z )=4 (1−x2 ) z−x

Y los valores iniciales:

t=0 , x=0.75 , z=0

Ahora procedemos a la ventana de comandos de Matlab, ingresando lo anteriormente mencionado:

>> f=@(t,x,z) z;

>> g=@(t,x,z) 4.*(1-x.^2).*z -x;

>> t=0;x=0.75;z=0;

>> h=0.1;m=50;

>> [u,v,w]=rk42(f, g, t, x, z, h, m);

Y de lo cual se obtiene el siguiente gráfico:

b) Realice alguna experimentación computacional cambiando el parámetro h=0.1, 0.05,0.025 y algún dato del problema. Escriba alguna observación acerca de la solución numérica calculada, tales como existencia, estabilidad, convergencia puntual, error de truncamiento.

Con h=0.1 y cambiando el valor incial de x=0.75 a 0.85Nos queda el mismo sistema pero varían las condiciones iniciales en la instrumentación computacional.

>> f=@(t,x,z) z;

>> g=@(t,x,z) 4.*(1-x.^2).*z -x;

>> t=0;z=0;

>> h=0.1;m=50;

>> x=0.75;[u,v,w]=rk42(f, g, t, x, z, h, m);x=0.85;[U,V,w]=rk42(f, g, t, x, z, h, m);

>> plot(u,v,'o','Color',[0,0,0],'MarkerEdgeColor','k','MarkerFaceColor','r','MarkerSize',5),hold on, grid on;

>> hold on, grid on,plot(U,V,'o','Color',[0,0,0],'MarkerEdgeColor','k','MarkerFaceColor','g','MarkerSize',5),hold on, grid on;

>>axis([ 0 1 0 1])

Los puntos rojos son de la solución original y los verdes son solución perturbada.

Ahora con h=0.05 y con los datos anteriores

>> f=@(t,x,z) z;>> g=@(t,x,z) 4.*(1-x.^2).*z -x;>> t=0;z=0;>> h=0.05;m=50;>> x=0.75;[u,v,w]=rk42(f, g, t, x, z, h, m);x=0.85;[U,V,w]=rk42(f, g, t, x, z, h, m);

>>plot(u,v,'o','Color',[0,0,0],'MarkerEdgeColor','k','MarkerFaceColor','r','MarkerSize',5),hold on, grid on;>> hold on, grid on,plot(U,V,'o','Color',[0,0,0],'MarkerEdgeColor','k','MarkerFaceColor','g','MarkerSize',5),hold on, grid on;>>axis([ 0 1 0 1])

Ahora con h=0.0025 y con los datos anteriores

>> f=@(t,x,z) z;>> g=@(t,x,z) 4.*(1-x.^2).*z -x;>> t=0;z=0;>> h=0.0025;m=1000;>> x=0.75;[u,v,w]=rk42(f, g, t, x, z, h, m);x=0.85;[U,V,w]=rk42(f, g, t, x, z, h, m);>>plot(u,v,'o','Color',[0,0,0],'MarkerEdgeColor','k','MarkerFaceColor','r','MarkerSize',5),hold on, grid on;>> hold on, grid on,plot(U,V,'o','Color',[0,0,0],'MarkerEdgeColor','k','MarkerFaceColor','g','MarkerSize',5),hold on, grid on;>>axis([ 0 1 0 1])

Conclusiones

Considerando el problema del tema 2, se puede decir que se trata de un problema correctamente planteado, ya que al realizar distintos cambios en las condiciones del problema y en la distancia entre los puntos utilizados por el método numérico, pudimos constatar de que la variación de la forma de la solución es la mínima cuando los cambios en sus condiciones son mínimos; es decir, que la solución no presenta puntos singulares conforme se varían sus parámetros, además tiende a cambiar lentamente, sin algún cambio brusco; mencionado lo siguiente, podemos asegurar que se trata de un método estable.Variando el valor de “h”, se pueden observar más puntos pertenecientes a la solución, lo cual nos ayudó a observar que la solución presenta convergencia. Aunque el parámetro perturbado fue el valor inicial de la función x, cuando perturbamos el valor del parámetro z, que representa la derivada de x con respecto a t, observamos que la solución tiende a cambiar bruscamente a la menor perturbación. Observando la gráfica 5, que es donde tenemos más puntos de las soluciones deseadas, notamos que el error de truncamiento de la solución perturbada se mantiene muy similar al error en la solución inicial.Finalmente, dado que tenemos convergencia en la solución, podemos afirmar que dicha solución existe.