Pràctica: INTEGRACIÓ SOBRE CORBES

• la integral d'una funció escalar f sobre una corba de parametritzada per

• i la circulació d'un camp vectorial al llarg d'una corba orientada,

i anàlogament per a corbes i camps vectorials a

En qualsevol d'aquests casos, hem d'acabar calculant una integral simple amb una funció que construïm a partir de les dades del problema. Per calcular aproximadament aquesta integral, podem usar diversos mètodes numèrics, alguns dels quals incorporats al Matlab. Com a mètode numèric molt simple, tenim la regla dels trapezis, donada per la funció trapz, però en general, obtindrem resultats millors aplicant algun dels mètodes de quadratura adaptativa com, per exemple, la quadratura adaptativa global que ve donada per la funció integral del Matlab (vegeu la descripció d'aquestes funcions a l'apèndix).

Camp vectorial i parametrització de la corba

Per il·lustrar el càlcul d'integrals sobre corbes, considerem el càlcul de l'àrea d'un domini a partir del teorema de Green, integrant el camp al llarg de la frontera del domini (recorreguda amb l'orientació compatible):

Com a exemple, calcularem l'àrea del domini definit en coordenades polars per és a dir, limitat per l' espiral logarítmica i el segment entre els seus punts inicial i final. Només caldrà integrar al llarg de l'espiral, ja que el segment és ortogonal al camp Compararem l'aproximació obtinguda amb el valor exacte:

Comencem definint el camp vectorial i la parametrització de la corba que recorre la frontera del domini en sentit antihorari (l'orientació compatible).

Dibuix de la corba i del camp vectorial

Comencem dibuixant la corba, usant plot, i tot seguit construïm una xarxa de punts en un rectangle de i dibuixem el camp vectorial a tots aquests punts amb la funció de Matlab quiver (a usaríem plot3 i quiver3).

close all %tanquem figures anteriors (optatiu)

figure, hold on, axis equal

%creem nova figura, sobre la qual podem anar afegint grafics,

%i igualem les proporcions dels eixos

 

N=100;

tv=linspace(0,2*pi,N+1); %rang de N+1 valors del parametre t

% Nota: alternativament podem generar el mateix rang de valors

% a partir de l'increment h del parametre de la corba,

% N=100;

% h=(2*pi-0)/N;

% tv=0:h:2*pi;

plot(x(tv),y(tv),'r') %dibuix de la corba

 

plot([x(0),x(2*pi)],[y(0),y(2*pi)],'r') %tambe dibuixem el segment

 

xi=linspace(-1.1*R,1.1*R,12);

yi=linspace(-1.1*R,1.1*R,12);

%per dibuixar el camp vectorial en una xarxa de 12x12 punts

%sobre el rectangle [-1.1*R,1.1*R]x[-1.1*R,1.1*R],

%comencem definint rangs de valors en les coordenades x,y

 

[xm,ym] = meshgrid( -1.1*R:R/5:1.1*R, -1.1*R:R/5:1.1*R );

%obtenim dues matrius xm i ym, que contenen la coordenada x

%i la coordenada y, respectivament, de tots els punts de la xarxa

 

Pm=P(xm,ym); %avaluem les components del camp vectorial a tots els punts

Qm=Q(xm,ym);

 

quiver(xm,ym,Pm,Qm) %podem afegir-hi un factor, "quiver(xm,ym,Pm,Qm,1.25)"

Càlcul de la integral usant la regla dels trapezis

sobre el rang de valors del paràmetre, obtenint un altre vector amb els valors Tot seguit apliquem la regla dels trapezis per tal d'obtenir una aproximació de la integral.

Càlcul de la integral usant quadratura adaptativa

En aquest cas hem de definir com una funció, a partir de les funcions P, Q, x, y, Dx, Dy ja definides.

Apèndix: Integració en 1 variable

Veurem com calcular una integral simple amb les funcions trapz i integral.

Escollim un nombre i dividim l'interval en N subintervals de longitud que anomenem el pas d'integració,

i calculem els valors Aproximem la integral de g sobre cadascun dels subintervals per l'àrea del trapezi limitat pel segment que uneix els punts i

Sumant aquests resultats per a obtenim una aproximació per a la integral que és la regla dels trapezis (o mètode dels trapezis):

Aquest mètode ens ve donat per la funció de Matlab trapz, que té com a arguments el vector o rang de valors format per les abscisses i un altre vector format per les ordenades

Com a exemple, calculem usant la regla dels trapezis amb subintervals. El valor exacte és

L'error en la regla dels trapezis és Això vol dir que si dividim el pas d'integració per 10, l'error es dividirà aproximadament per 100 (és a dir, guanyem 2 xifres decimals correctes). Podem verificar aquest fet a l'exemple anterior, repetint els càlculs amb subintervals.

Diverses funcions de Matlab incorporen mètodes d'integració més acurats que la regla dels trapezis, que podem aplicar si coneixem una fórmula explícita de la funció a integrar. En particular, podem usar integral (quadratura adaptativa global), en què no s'ha d'especificar un pas d'integració (que serà variable) sinó toleràncies per a l'error absolut i relatiu. Així, escriurem q=integral(g,a,b, 'AbsTol',tol1, 'RelTol',tol2) per a calcular una integral amb tolerància tol1 per a l'error absolut i tolerància tol2 per a l'error relatiu (cal dir però que això és orientatiu i depèn en bona part de les propietats de la funció que volem integrar). Si no s'especifiquen les toleràncies, per defecte són 1e-10 i 1e-6 respectivament.

Il·lustrem-ho amb la integral primer amb les toleràncies per defecte i després escollint toleràncies tol1=tol2=1e-12. El valor exacte és

Nota: Una opció alternativa per definir la funció a integrar és, en comptes d'usar una ordre del tipus g=@(t)..., fer-ho a través d'un fitxer, anomenat per exemple gfun.m, que contingui les línies següents,