Matlab - Reducir tiempo de ejecución de las matrices

   
Vista:
Imágen de perfil de Alex

Reducir tiempo de ejecución de las matrices

Publicado por Alex (1 intervención) el 05/06/2016 23:32:13
Estaba realizando un programa de diferencias finitas, pero cada vez que incremento el valor de la variable n, el programa se demora demasiado en mostrar el resultado, lo cual no es conveniente ya que necesito que el programa se ejecute en el menor tiempo posible.


1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
a=0;%limite inferior del intervalo
b=2;%limite superior del intervalo
n=10000;%numero de subdivisiones
h=b/n;%incremento
syms x;
%Aplicaremos el método de diferencias dinitas a la ecuacion diferencial 
%-d/dx(p du/dx)+qu=f, donde u(0)=0=u(5)
%Donde p(x)=1+x, q(x)=x
f=x*(-1+x*(-2+x))-3*x;%f(x)=40*x+10*x^2 *(2-x)
p=1+x;%p(x)=1+x
q=x;%q(x)=x
ua=0;% u(0)=0
ub=0;%u(b)=0
%Se formará el sistema de ecuaciones Au=b
% A = zeros(n-1);%Matriz tridiagonal, donde la diagonal principal tiene 
%elementos de la forma (p(xi^)+ p(x_i+1^))/h^2 + q(xi), el las diagonales superior e inferior
%tienen la forma - p(x_i+1^)/h^2, con i=1,...,n-1
x=a:h:b;%estos son los xj=j*h
xjt=zeros(n,1);%Creamos un vector que contengan los xj^=(x_j-1 + xj)/2
               % y x_j+1^ = (xj + x_j+1)/2
for i=2:length(xjt)
   xjt(i-1)=(x(i-1)+x(i))/2; %xj^=(x_j-1 + xj)/2
   xjt(i)=(x(i)+x(i+1))/2;%  x_j+1^ = (xj + x_j+1)/2
end
 
p1=inline(p);%Sea evalua la funcion p en  xj^
p1(xjt);%Se calcula los p(xj^), con j=1,...n-1;
q1=inline(q) ;% Sea evalua la funcion q en los puntos xj
q1(x);%Se calcula los q(xj), con j=1,n
%%%%%%%%%%%%%%%%%%%
bt= zeros (n-1,1);%Es la matriz b del sistema de ecuaciones Au=b%
sm=eval(f);%Evaluamos la función f en los puntos xj
for i=1:n-1
    bt(i)=sm(i); %Los elementos son de la forma f(xj), j=1,...,n-1
end
B=zeros(n-1,3);%Arreglo de nx3, donde solo se han puesto los valores 
%no nulos de la matriz A
for i=1:n-1
    for j=1:3
        if ((i==1) && (j==1))||((i==n) && (j==3))
            B(i,j)=0;
        end
        if (i>= 2)
            B(i,1)= -p1(xjt(i))/h^2;
            B(i,3)=-p1(xjt(i+1))/h^2;
        end
        B(i,2)= (p1(xjt(i))+p1(xjt(i+1)))/h^2 + q1(x(i+1));
    end
end
% Como el arreglo B es de nx3, el ultimo elemento del vector sol no se puede
% obtener por la eliminacion gaussiana aplicada directamente en el arreglo B
for i=2:length(bt)
temp=B(i-1,3)/B(i-1,2);
B(i,2)=B(i,2)-temp*B(i,1);
B(i,1)=0;
bt(i)=bt(i)-temp*bt(i-1);
end
%Sustitución hacia atrás
sol=zeros(n-1,1);
sol(n-1)=bt(n-1)/B(n-1,2);%El ultimo elemento de bt 
for j=length(bt)-1:-1:1
    sol(j)=(bt(j)-B(j,3)*sol(j+1))/B(j,2);
end
plot (sol)
Valora esta pregunta
Me gusta: Está pregunta es útil y esta claraNo me gusta: Está pregunta no esta clara o no es útil
0
Responder

Reducir tiempo de ejecución de las matrices

Publicado por crs (13 intervenciones) el 07/06/2016 17:37:44
Hola,

Lo que más suele consumir en este tipo de códigos son los bucles for, especialmente los anidados. Si quieres reducir el tiempo, deberías intentar sustituirlos utilizando el formato matricial para acelerar los cálculos. Por ejemplo, el bucle de las líneas 36 a 49 puedes sustituirlo por el siguiente código, que debería obtener el mismo resultado para B, pero en un tiempo mucho menor:

1
2
3
4
5
6
B = zeros(n-1,3);
B(1,1) = 0;
B(1,2) = (p1(xjt(1))+p1(xjt(2)))/h^2 + q1(x(2));
B(2:end,1) = -p1(xjt(2:n-1))/h^2;
B(2:end,2) = (p1(xjt(2:n-1))+p1(xjt(3:n)))/h^2 + q1(x(3:n)');
B(2:end,3) = -p1(xjt(3:n))/h^2;

Un saludo,
Valora esta respuesta
Me gusta: Está respuesta es útil y esta claraNo me gusta: Está respuesta no esta clara o no es útil
1
Comentar