Introducción informal a Matlab y Octave - Ejercicios resueltos (I)

Metodología de programación

8.1  Ejercicio. Cálculo de un gradiente

8.1.1  Guía para la resolución del ejercicio

8.1.2  Solución del ejercicio

function[grad]=mygradient(M)
 
     grad.x = 0.5*(M(:,3:end)-M(:,1:end-2));
 
     grad.y = 0.5*(M(3:end,:)-M(1:end-2,:));
  

8.2  Ejercicio. Diseño de una tobera

1. Resolver el sistema de ecuaciones anterior y representar las soluciones de P=0.9, P=1 y P=1.2

8.2.1  Guía para la resolución del ejercicio

octave:1> fplot('[x. ^2+0.1,sqrt(x)]',[0,1])
  

octave:2> help fsolve
  

8.2.2  Solución del ejercicio

function[f]=ttb(x)
     global Pi
     f(1)=x(1)/(2*sqrt(x(1)*x(3)+x(2)))-2*Pi*x(3);
     f(2)=Pi*x(3)**2+0.1-sqrt(x(1)*x(3)+x(2));
     f(3)=sqrt(x(1)+x(2))-1;
 end   
   
 function[f]=tobera(x,a,b,l,P)
     if x<l
       f=P*(x*x)+0.1;
     else
       f=sqrt(a*x+b);
     endif
 end   
 x0=[1 1 1];
 P=[0.9 1 1.2];
 hold on   
 for i=1:3
   global Pi
   Pi=P(i);
   [res]=fsolve('ttb',x0);
   xcoord=linspace(0,1,100);   
   for j=1:100
     tob(j)=tobera(xcoord(j),res(1),res(2),res(3),P(i));
   end
   plot(xcoord,tob)
 end   
    
 hold off
  

---

Figure 8.1: Resultado del script

---

8.3  Ejercicio. El atractor de Lorentz

8.3.1  Guía para la resolución del Ejercicio

• Primero escribimos la función correspondiente a la ecuación diferencial. En las subrutinas de resolución de EDOs siempre tenemos que introducir la función de la forma dx/dt=f(x,t) pudiendo ser cualquiera de las variables un vector.
• La rutina de Octave que nos integra la ecuación diferencial es lsode; en Matlab utilizaremos ode45 porque el problema no es stiff. Utilizar la ayuda para saber de qué manera tenemos que escribir la función de la ecuación diferencial.
• El comando para representar curvas paramétricas en 3 dimensiones es plot3.
• Escribiremos todo en un script y lo ejecutaremos

8.3.2  Solución del ejercicio

8.3.2.1  Octave

x=0;
 function xdot=func(x,t)
 a=10;b=28;c=8/3;
 xdot(1,1)=a*(x(2)-x(1));
 xdot(2,1)=x(1)*(b-x(3))-x(2);
 xdot(3,1)=x(1)*x(2)-c*x(3);
 end
 
 x0=[1;1;1];
 t=linspace(0,50,5000);
 tic;x=lsode( "func ",x0,t);toc
 plot3(x(:,1),x(:,2),x(:,3))
 

---

Figure 8.2: Resultado del script

---

8.3.2.2  Octave no stiff

lsode_options('integration method','non-stiff')
  

8.3.2.3  Octave y C++

t=linspace(0,50,5000);
 tic;x=lsode( "eqlorentz ",[1;1;1],t);toc
 plot3(x(:,1),x(:,2),x(:,3))
 

8.3.2.4  Octave y C++ no stiff

8.3.2.5  Octave, C++ y Fortran

8.3.2.6  Matlab

function xdot=lorentz(t,x)
   a=10;b=28;c=8/3;
   xdot(1,1)=a*(x(2)-x(1));
   xdot(2,1)=x(1)*(b-x(3))-x(2);
   xdot(3,1)=x(1)*x(2)-c*x(3);
 end
  

x=0;t=0;
 tic;[t,x]=ode45(@lorentz,[0,50],[1,1,1]);toc
 plot3(x(:,1),x(:,2),x(:,3))
  

8.4  Ejercicio. Cálculo de una integral doble

8.4.1  Guía para la resolución del ejercicio

 >> help quad2dg
 quad2dg is the user-defined function from the file
 /usr/share/octave/2.1.64/site/m/octave-forge/integration/quad2dg.m   
 usage:  int = quad2dg('Fun',xlow,xhigh,ylow,yhigh)
 or
        int = quad2dg('Fun',xlow,xhigh,ylow,yhigh,tol)   
 This function is similar to QUAD or QUAD8 for 2-dimensional integration,
 but it uses a Gaussian quadrature integration scheme.
       int     -{}- value of the integral
       Fun     -{}- Fun(x,y) (function to be integrated)
       xlow    -{}- lower x limit of integration
       xhigh   -{}- upper x limit of integration
       ylow    -{}- lower y limit of integration
       yhigh   -{}- upper y limit of integration
       tol     -{}- tolerance parameter (optional)   
   
 Note that if there are discontinuities the region of integration
 should be broken up into separate pieces.  And if there are singularities,
 a more appropriate integration quadrature should be used
 (such as the Gauss-Chebyshev for a specific type of singularity).
  

8.4.2  Solución del ejercicio

• Con inline
  • Matlab:
    
    

    >> dblquad(inline('exp(-(x.*x+y.*y))'),-10,10,-10,10)
     
     ans = 3.1416
     
     

  • Octave:
    
    
    
    

    >> quad2dg(inline('exp(-(x.*x+y.*y))'),-10,10,-10,10)
     
     ans = 3.1416
     
     

• Con una función anónima:
  • Matlab:
    
    
    
    

    >> dblquad(@(x,y) exp(-(x. ^2+y. ^2)),-10,10,-10,10)
     
     ans = 3.1416
     
      

    
    
    
  • Octave:
    
    
    
    

    >> quad2dg(@(x,y) exp(-(x. ^2+y. ^2)),-10,10,-10,10)
     
     ans = 3.1416