Derivaci n e Integraci n num rica o o e by rockman16

VIEWS: 0 PAGES: 5

									               o             o     e
       Derivaci´n e Integraci´n num´rica
                                           a
                      Complementos de Matem´tica Aplicada
                                    Curso 2007-2008


1.             o     e
       Derivaci´n num´rica
1.1.   Polinomios
   Para derivar un polinomio p(x) = a1 xn + ... + an+1 , dado en un vector p, usaremos
q=polyder(p), que nos dar´ el polinomio q(x) = na1 xn−1 + ... + an .
                         a

Ejercicio 1 Encuentra y clasifica los extremos relativos de p(x) = x4 − x2 + 1. Para ello:
  1. Encuentra q(x) = p (x).
                       ı               a
  2. Halla todas las ra´ces de q(x). Gu´rdalas en un vector r.
                          o
  3. Halla s(x) = p (x) (l´gicamente s(x) = q (x)).
         u
  4. Eval´a el polinomio s(x) en los valores guardados en r.

1.2.   Splines
                                                o
   Para derivar un spline obtenido mediante las ´rdenes csape o spline, se utiliza el comando
fnder.
Ejemplo 1 El siguiente conjunto de instrucciones genera un spline pp que interpola un
conjunto de puntos aleatorio y su derivada, y pinta ambas funciones:
% Generamos una tabla de 11 puntos aleatoria
x=0:10;
y=rand(size(x));

% Generamos el spline
pp=csape(x,y);

% Lo pintamos
xx=0:.01:10;           % Para ello tomamos 1001 puntos en [0,10]
yy=ppval(pp,xx);       % Evaluamos el spline en esos 1001 puntos
plot(xx,yy)
hold

% Generamos su derivada
ppd=fnder(pp);

% La pintamos
yyd=ppval(ppd,xx);     % Para ello evaluamos la derivada en los 1001 puntos de antes
plot(xx,yyd,’r’)       % Pintamos en rojo: ’r’

                                              1
                          a             ı
Ejercicio 2 Halla los m´ximos y los m´nimos relativos del spline generado en el apartado 1.
                                                                  ı
del ejercicio 3 del tema Datos y Polinomios. Para encontrar las ra´ces de un polinomio a trozos,
usa fnzeros. Mira en la ayuda doc fnzeros


1.3.    Otras funciones
                             o                           e                                 o
    Matlab no trae una funci´n general para evaluar num´ricamente la derivada de una funci´n
                               o                                             o
en un punto. La siguiente funci´n calcula la derivada de orden o de una funci´n dada en f en
                                  o
un punto x. Para ello calcula la f´rmula de diferencias centradas de orden 2n con paso h. El
  e
m´todo consiste en calcular el polinomio interpolador en los nodos correspondientes, derivarlo
                            e
y evaluar en el punto x. Obs´rvese el uso de argumentos de entrada opcionales.
    http://enol.epsig.uniovi.es/complementos/index_archivos/der.m

function y=der(f,x,o,n,h)
            o      e
% Derivaci´n num´rica
%
%      a         o
    Est´ funci´n que calcula la derivada de orden o de
%               o
    una funci´n dada en f en un punto x. Para ello
%                  o
    calcula la f´rmula de diferencias centradas de orden
%   2n con paso h.
%   Uso
%
%      y=der(f,x<,o,n,h>)
%
%   Los argumentos o, n y h son opcionales, siendo
%   sus valores por defecto o=1, n=2 y h=0.001.
%
%                         o
    Ejemplos de utilizaci´n
%
%        d1=der(@f1,0)
%        d2=der(@f1,0,2)
%        d3=der(@f1,0,2,2)
%        d4=der(@f1,0,2,2,1e-4)
%
%             o
    La funci´n ha de poder evaluarse sobre un vector:
%                            o
    Correcto: x.^6 (evaluaci´n sobre un vector)
%                              o
    Incorrecto: x^6 (evaluaci´n sobre un punto)
%
%         e                           o
    Tambi´n pueden usarse funciones an´nimas
%
%        d5=der(@(x) sin(x),0)
%
if nargin<=4
    h=1e-3;
end
if nargin<=3
    n=1;

                                               2
end
if nargin==2
    o=1;
end
X=x-n*h:h:x+n*h;
Y=f(X);
                                        a
[p,S,mu]=polyfit(X,Y,2*n); %Una manera m´s complicada de usar polyfit
                             n                   e
                           %A~ade estabilidad num´rica, pero ahora p
                           %ya no es un polinomio en la variable x, sino
                           %en la variable (x-mu(1))/mu(2)
for k=1:o
p=polyder(p)/mu(2);        %Hay que tener cuidado al derivar por el cambio
                           % de variable
end
y=polyval(p,(x-mu(1))/mu(2)); %Y lo mismo al evaluar el polinomio.



Ejercicio 3 Calcular a mano las derivadas de las siguientes funciones en los puntos que
                                 o                                           e
indican. Utilizar luego la funci´n der para calcular las mismas derivadas num´ricamente.
                                       a
Probar con distintos valores de los par´metros n y h en cada caso.

                               f (x) = x6 + 3x2 + 1, en x = 0

                                g(x) = sin(300x), en x = 2π

                            k(x) = |x|, en x = 0 y en x = 10−6


2.              o     e
       Integraci´n num´rica
2.1.    Polinomios
   Para hallar una primitiva del polinomio p(x) = a1 xn + ... + an+1 , dado en un vector p,
                                                         a1
usaremos q=polyint(p), que nos dar´ el polinomio q(x) = n+1 xn+1 + ... + an+1 x.
                                   a


                                                  1
Ejercicio 4 Sea p(x) = x4 − x2 + 1. Calcula 0 p(x)dx. Para ello calcula una primitiva q(x)
             u
de p(x), eval´ala en los extremos del intervalo y aplica la regla de Barrow.


2.2.    Splines
                                                                    o
    Para hallar una primitiva de un spline pp obtenido mediante las ´rdenes csape o spline,
se utiliza el comando ppi= fnint(pp).

                                                                       6,1
Ejercicio 5 Sea pp(x) el spline obtenido en el ejercicio 2. Calcula   2,5
                                                                             pp(x)dx.

                                              3
2.3.    Otras funciones
                                                                e
    Matlab 7 trae varias funciones para hacer integrales num´ricamente. Todas subdividen el
                                                              o
intervalo tantas veces como sea necesario y realizan integraci´n adaptativa. Las principales son
                                           o
quad, que usa Simpson y quadl, que usa f´rmulas adapativas de Gauss-Lobatto.
                                   e      a                                  o
    A partir de Matlab 7.5, tambi´n est´ disponible quadgk, que usa la f´rmula de Gauss–
Kronrod (7-15).
                                                                    o
    Para usar ambas funciones es necesario proporcionar una funci´n que pueda ser evaluada
                                   e
en un vector de nodos. El uso es id´ntico en ambos casos:

>> I=quad(@f,a,b)

>> I=quadl(@f,a,b)}
             b
hacen I =   a
                               e                             o
                 f (x)dx. Tambi´n se pueden usar funciones an´nimas.
                             1 −x2
Ejemplo 2 Para hacer        0
                              e dx,                    ı
                                           la orden ser´a

>> I=quadl(@(x) exp(-x.^2),0,1)

   Observa que si no pones el punto antes de la potencia, da un error.

                               o           e                        e
Mediante el uso de funciones an´nimas tambi´n se pueden integrar num´ricamente polinomios
y splines:

>> I=quadl(@(x) polyval(p,x),0,1)

                                                     e
>> I=quad(@(x) ppval(pp,x),0,10) % ¡Es exacta! ¿Porqu´?

>> I=quad(@ppval,0,10,[],[],pp) % Otra manera debido a la estructura especial
                              % de los splines

                                  e
Ejercicio 6 Realizar por los dos m´todos (Simpson y Gauss-Lobato)
                            1                    2π             1                 1
                                f (x)dx,              g(x)dx,        k(x)dx           x10 dx
                           −1                0                  −1            0

¿En cuales se puede asegurar que se ha obtenido el resultado exacto?



                          o                                 o
Ejercicio 7 En la versi´n 7.0 de Matlab, la programaci´n de quad y quadl tiene un peque˜o      n
                                                                       o
truco para tratar los puntos singulares. Si se encuentra con una divisi´n por 0 la evita moviendo
un poco el punto. Esto se hace de forma transparente para el usuario, simplemente da un aviso
y calcula la integral. Calcula
                                        1              1
                                           1             1
                                          √ dx, y          dx.
                                      0     x         0 x

                                                                                    e
Observa que en la primera sale un aviso y en la segunda dos avisos distintos. ¿Porqu´?



                                                           4
                  a                            o
Ejercicio 8 En mec´nica es frecuente la aparici´n de integrales del tipo
                                  α
                                           dx
             F (α, k) =                                               ı
                                                        , (Integral el´ptica de primera clase)
                              0        1 − k 2 sin2 x
y
                              α
            E(α, k) =                 1 − k 2 sin2 xdx, (Integral el´ptica de segunda clase)
                                                                    ı
                          0

                       o                                                           o
    a) Crear una funci´n function y=F(alfa,k) que calcule F (α, k), y una funci´n function
                                                               o
y=E(alfa,k) que calcule E(α, k). Usar para ello quadl. N´tese que en cada caso hay que
                   o                           a
integrar una funci´n que depende de un par´metro (k). O bien se crea otro fichero de funci´n  o
con el integrando y se pasa k como variable global o bien se escribe el integrando como una
      o    o                                 a
funci´n an´nima (lo cual resulta mucho m´s sencillo).
               e                                                               a
    b) En un p´ndulo simple de longitud , bajo gravedad g, lanzado desde un ´ngulo θ desde la
                                                                  o                     o
vertical con velocidad angular inicial ω, el periodo P de oscilaci´n viene dado por la f´rmula
                                                        4   π
                                                P =       F   ,k ,
                                                        β   2
donde
                                                                       ω2
                                                   1 − cos θ −         2β 2
                                           k=
                                                               2
y
                                                               g
                                                    β=             .

   Tomando g = 9,8m/s2 , = 1m, θ = π/2 y ω = 0, calc´lese P .
                                                           u
                       o
   c) Crea una funci´n function P=periodo(theta) que nos proporcione el periodo en
     o      a                           a
funci´n del ´ngulo inicial. Dibuja su gr´fica en [0, π/2]. Como seguramente no has programado
periodo.m de manera vectorial, usa el siguiente programa:

th=linspace(0,pi/2,200);                  % Definimos las amplitudes
P=zeros(size(th));                        % Inicializamos el vector de periodos
for j=1:length(th)
    P(j)=periodo(th(j));                  % Calculamos los periodos uno a uno
end
plot( th*180/pi, P)                       % Dibujamos grados vs periodos

Haz un zoom para ver la gr´fica entre 0 y 10o .
                          a
               e                                          a
    d) Busca “p´ndulo” en la wikipedia y lee el segundo p´rrafo. ¿Es eso cierto tal y como lo
       ı?          a                        o           ı             u
dice ah´ Mira la gr´fica del apartado c). ¿C´mo lo habr´as redactado t´?




                                                           5

								
To top