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
0
