TRABAJO FINAL DIMA 4° 01 JUAN ANTONI ANTONIO O GONZALEZ GONZALEZ LOPEZ
Método de Cuadratura de Gauss Gauss En análisis análisis numérico una cuadratura es una aproximación aproximación de una integral definida de una función. Una cuadratura de Gauss n-puntos n-puntos llamada llamada así por Carl Friedrich Friedrich Gauss. Gauss. Es una cuadratura que selecciona los puntos de la evaluación de manera óptima y no en una forma igualmente espaciada, construida para dar el resultado de una polinomio de grado 2n-1 o menos, elegibles para los puntos x i y los coeficientes w i para i=1,...,n. El dominio de tal cuadratura por regla es de [−1, 1] dada por :
Tal cuadratura dará resultados precisos solo si f(x) es aproximado por un polinomio dentro del rango [−1, 1]. Si la función puede ser escrita como f(x)=W(x)g(x), donde g(x) es un polinomio aproximado y W(x) es conocido.
Programa function y = gaussquadratura( f, a, b, TOL ) h2 = ( b - a ) / 2; fm = feval ( f, a + h2 - sqrt(1/3)*h2 ); fp = feval ( f, a + h2 + sqrt(1/3)*h2 ); sab = h2 * ( fm + fp ); [approx eest nfunc] = as ( sab, f, a, b, TOL ); if (( nargout == 0 ) if s = sprintf ( '\t\t approximate value of integral: \t %.12f \n', \n' , approx ); s = sprintf ( '%s \t\t error estimate: \t\t\t\t\t %.4e \n' , s, eest ); s = sprintf ( '%s \t\t number of function evaluations: \t %d \n' , s, nfunc+2 ); disp ( s ) else y = approx; end;; end return function [app, est, nf] = as ( sab, f, a, b, TOL ) h2 = ( b - a ) / 4; fml = feval ( f, a + h2 - sqrt(1/3)*h2 ); fpl = feval ( f, a + h2 + sqrt(1/3)*h2 ); fmr = feval ( f, b - h2 - sqrt(1/3)*h2 ); METODOS NUMERICOS
PROF: ESTEBAN ISIDRO PIOQUINTO Página 1
fpr = feval ( f, b - h2 + sqrt(1/3)*h2 ); sac = h2 * ( fml + fpl ); scb = h2 * ( fmr + fpr ); errest = abs ( sab - sac - scb );
if ( errest < (10.0*TOL) ) app = sac+scb; est = errest / 10.0; nf = 4; return; else [a1 e1 n1] = as ( sac, f, a, a+2*h2, TOL/2.0 ); [a2 e2 n2] = as ( scb, f, a+2*h2, b, TOL/2.0 ); app = a1 + a2; est = e1 + e2; nf = n1 + n2 + 4; return; end; Ejemplo en matlab >> a=2 %es el valor inferior de la integracion a= 2 >> b=3 % es el valor mayor de la integracion b= 3 >> TOL=0.1 % es la tolerancia del error de convergencia TOL = 0.1000 >> f=@(x)(1.44./x.^2)+0.24 % es el nombre de la integral f= @(x)(1.44./x.^2)+0.24 >> y = gaussquadrature2( f, a, b, TOL ) y=
0.4800
MÉTODO TRAPEZOIDAL La regla del trapecio o regla trapezoidal es la primera de las fórmulas cerradas de Newton-Cotes. Corresponde al caso en donde el polinomio de aproximación es de primer orden.
En donde f1(x) corresponde a una línea recta que se representa como:
El área bajo la línea recta es una aproximación de la integral de f(x) entre los límites a y b:
El resultado de la integración es:
Programa function y = trapezio(valor min,valor max, interv,fun ) y = 0; step = (valor max - valor min) / interv; for j = valor min : step : valor max y = y + feval(fun,j); end; y = (y - (feval(fun, valor min) + feval(fun, valor max))/2) * step;
Programa de la función
function y = fun_for_integration(x) y = x^3; ejemplo en matlab >> z1 = trapezio(0, 2, 10, 'fun_for_integration') z1 = 4.0400 >> z2 = trapezio(0, 3, 65, 'fun_for_integration') z2 = 20.2548
Método de Simpson En análisis numérico, la regla o método de Simpson (nombrada así en honor de Thomas Simpson) y a veces llamada regla de Kepler es un método de integración numérica que se utiliza para obtener la aproximación de la integral:
.
Derivación de la regla de Simpson Consideramos el polinomio interpolante de orden dos P2(x), que aproxima a la función integrando f(x) entre los nodos x 0 = a, x1 = b y m = (a+b)/2. La expresión de ese polinomio interpolante, expresado a través de la Interpolación polinómica de Lagrange es:
Así, la integral buscada se puede aproximar como:
Programa function Int = Simpson(x,y) foo = 1:length(x); fin = length(x); h = (x(fin)-x(1))/(fin-1); odd = find(rem(foo,2)~=0);
even = find(rem(foo,2)==0); IntOdd = 2*y(odd); IntEven = 4*y(even); if rem(foo(fin),2)~=0 Int = (h/3)*(sum(IntOdd) + sum(IntEven)-y(1)-y(fin)); else Int = (h/3)*(sum(IntOdd) + sum(IntEven)-y(1)-3*y(fin)); end end Ejemplo en matlab >> x=1:5 x= 1
2
3
4
5
7
8
9
>> y=5:10 y= 5
6
10
>> Int = Simpson(x,y) Int = 28
Método de Euler La fórmula o relación de Euler, atribuida a Leonhard Euler, establece que:
para todo número real x . Aquí, e es la base del logaritmo natural, i es la unidad imaginaria y senx y cosx son funciones trigonométricas. Ó bien:
siendo Z la variable compleja formada por : Z=a+ix.
Demostración usando las Series de Taylor La fórmula de Euler ilustrada en el plano complejo. Sabiendo que:
y así sucesivamente. Además de esto, las funciones e x , cos( x ) y sin( x ) (asumiendo que x sea un número real) pueden ser expresadas utilizando sus series de Taylor alrededor de cero.
Programa
function [t,y]=feuler(f,t0,b,y0,n) h=(b-t0)/n; t=zeros(1,n+1); y=zeros(1,n+1); index=[0:1:n]; t=t0+h*index; y=y0; for i=2:n+1 y(i)=y(i-1)+h*feval(f,t(i-1),y(i-1)); end Ejemplo en matlab
>> f=@(y,t)5*y./(t+1)-y f= @(y,t)5*y./(t+1)-y >> [tvals,yvals]=feuler(f,0,4,1,20) tvals = Columns 1 through 9
0 0.2000 0.4000 0.6000 0.8000 1.0000 1.2000 1.4000 1.6000 Columns 10 through 18 1.8000 2.0000 2.2000 2.4000 2.6000 2.8000 3.0000 3.2000 3.4000 Columns 19 through 21 3.6000 3.8000 4.0000 yvals = Columns 1 through 9 1.0000 1.0000 1.0600 1.1742 1.3301 1.5135 1.7113 1.9139 2.1144 Columns 10 through 18 2.3081 2.4922 2.6649 2.8252 2.9726 3.1071 3.2289 3.3383 3.4359 Columns 19 through 21 3.5224 3.5984 3.6648 >> [t,y]=feuler(f,0,4,1,20); >> yy=(t+1).^5.*exp(-t); >> plot(t,y,'x',t,yy)
0
0
0
0
0
Metodo de Runge-Kutta 2° Orden
0.5
1
1.5
2
2.5
3
3.5
La versión de segundo orden de la ecuación y i + 1 = y i + φ( x i, y i ,h)h es y i + 1 = y i + (a1k 1 + a2k 2)h
donde k 1 = f ( x i, y i) k 2 = f ( x i + pi h,y i + q11k 1h)
Los valores de a1, a2, p1 y q11 son evaluados al igualar el término de segundo orden de la ecuación dada con la expansión de la serie de Taylor. Se desarrollan tres ecuaciones para evaluar cuatro constantes desconocidas. Las tres ecuaciones son: a1 + a2 = 1 a2 p1 = 1 / 2 a2q11 = 1 / 2
Como se tienen tres ecuaciones con cuatro incógnitas se tiene que suponer el valor de una de ellas. Suponiendo que se especificó un valor para a2, se puede resolver de manera simultánea el sistema de ecuaciones obtenido: a1 = 1 − a2
Como se puede elegir un número infinito de valores para a2, hay un número infinito de métodos Runge-Kutta de segundo orden. Cada versión podría dar exactamente los mismos resultados si la solución de la EDO fuera cuadrática, lineal o una constante. Programa function[y]= rungekutta2(h) t=10; y(1)=0; n=10/h; a=0.2; b=0.4; x=0:h:t; for i=2:n+1; k1=h*(a-b*y(i-1)); k2=h*(a-b*(y(i-1)+k1)); y(i)=y(i-1)+0.5*(k1+k2); end; plot(x,y,'o') hold on Ejemplo en Matlab
>> h=0.5 h= 0.5000 >> [y]= rungekutta2(h) y= Columns 1 through 9 0 0.0900 0.1638 0.2243 0.2739 0.3146 0.3480 0.3754 0.3978 Columns 10 through 18 0.4162 0.4313 0.4436 0.4538 0.4621 0.4689 0.4745 0.4791 0.4829 Columns 19 through 21 .
0.4860 0.4885 0.4906
.45
.
.35
.3
.
.2
.15
.1
.05
0
1
2
3
4
5
6
7
8
9
El método de Runge-Kutta de cuarto orden Un miembro de la familia de los métodos Runge-Kutta es usado tan comúnmente que a menudo es referenciado como “RK4” o como “el método Runge-Kutta”. Definamos un problema de valor inicial como:
Entonces el método RK4 para este problema está dado por la siguiente ecuación:
Donde
10
Así, el siguiente valor (yn+1) es determinado por el presente valor (yn) mas el producto del tamaño del intervalo (h) por una pendiente estimada. La pendiente es un promedio ponderado de pendientes: k 1 es la pendiente al principio del intervalo; k 2 es la pendiente en el punto medio del intervalo, usando k 1 para determinar el valor de
y en el punto
usando el método de Euler
k 3 es otra vez la pendiente del punto medio, pero ahora usando k 2 para determinar el valor de y k 4 es la pendiente al final del intervalo, con el valor de y determinado por k 3
Promediando las cuatro pendientes, se le asigna mayor peso a las pendientes en el punto medio:
Esta forma del método de Runge-Kutta, es un método de cuarto orden lo cual significa que el error por paso es del orden de O(h5), mientras que el error total acumulado tiene el orden O(h4).
Programa function[y,z]=rungekutta4(h,t,y,z,r,a,b) n=r/h; x=0:h:t; for i=1:n; k1=h*(z(i)); l1=h*(-a*z(i)-b*y(i)); k2=h*(z(i)+l1/2); l2=h*(-a*(z(i)+l1)-b*(y(i)+k1/2)); k3=h*(z(i)+l2/2); l3=h*(-a*(z(i)+l2)-b*(y(i)+k2/2));
k4=h*(z(i)+l3/2); l4=h*(-a*(z(i)+l3)-b*(y(i)+k3/2)); y(i+1)=y(i)+((1/6)*(k1+(2*k2)+(2*k3)+k4)); z(i+1)=z(i)+((1/6)*(l1+(2*l2)+(2*l3)+l4)); end; hold on plot(x,y,'-') plot(x,z,'r-') Ejemplo en Matlab
.8
.6
.
.2
- .2
.4
.6
- .8 0
5
10
15
20
25
30
35
40
45
50
