Matlab - Matlab lento

% Discretizaci?n del problema Dx(a(x,t)DxU(x,t))=DtU(x,t)-f(x,t)

%function [A,f,ateorica,uteorica,fvector,avector,x,t]=calorparabolica

function [A,ateorica,uteorica,fvector,avector,x,t]=calorparabolica

clc; clear

J=64;N=64;

%xn=linspace(0,1,J);

%tn=linspace(0,1,N);

%parpool

xn=linspace(0,1,J*N);

tn=linspace(0,1,N*J);

syms xx tt

us = exp(tt).*sin(pi.*xx);

ps = 40.*xx.^2./(40.*xx.^2+1)+1+0.*tt;

ffs=-diff(ps.*diff(us,xx),xx)+diff(us,tt);

%ezplot(ps,[0,1])

%[xm,tm]=meshgrid(xn,tn);

%xm=xm';

%tm=tm';

x=xn;

t=tn;

%numericos

fn=subs(ffs,{xx,tt},{xn,tn});

un=subs(us,{xx,tt},{xn,tn});

pn = subs(ps,{xx,tt},{xn,tn});

%f=fn;

%plot(pn)

A = zeros(J*N,J*N);

%A = zeros(J+N,J+N);

% pasar a avector las matrices

uu =un(:);

fvector = fn(:);

%avector =pn(:);

avector = pn;

h=1/(J-1);

ateorica = pn;

uteorica = un;

%% Matriz

 for i = 1:J*N-2

%for i = 1:J+N

   if i == 1 % t < h/2 (0,h/2)

       A(i,i:i+1) = [((uu(i+1) - uu(i+2))) (5*uu(i+1) - 2*uu(i+2) - 3*uu(i))]/(h^2);

   elseif i == J+N % t > 1-h/2 (1-h/2,1)

       A(i,i-1:i) = [(5*uu(i-1) - 2*uu(i-2) - 3*uu(i))  (uu(i-1)-uu(i-2))]/(h^2);

   else % h/2 < t < 1-h/2 (h/2,1-h/2)

       A(i,i-1:i+1) = [(-uu(i) + uu(i-1)) (uu(i+1)-2*uu(i)+un(i-1)) (uu(i+1)-uu(i))]/(2*h^2);

   end

 end

 aa = pinv(A)*fvector;

 %delete(gcp('nocreate'))

 plot(aa)

% Discretizaci?n del problema Dx(a(x,t)DxU(x,t))=DtU(x,t)-f(x,t)

%function [A,f,ateorica,uteorica,fvector,avector,x,t]=calorparabolica

function  [A,ateorica,uteorica,fvector,avector,x,t]=calorparabolica

close all

J=30;

%funcion a la perfeccion hasta J=30, pero puedes intentar con otro

%números mayores

N=J;

%xn=linspace(0,1,J);

%tn=linspace(0,1,N);

%parpool

xn=linspace(0,1,N*J);

tn=linspace(0,1,J*N);

syms xx tt

us = vpa(exp(tt).*sin(pi.*xx),3);

ps = vpa(40*xx.^2./(40*xx.^2+1)+1+0.*tt,3);

ffs=vpa(-diff(ps.*diff(us,xx),xx)+diff(us,tt),3);

%ezplot(ps,[0,1])

%[xm,tm]=meshgrid(xn,tn);

%xm=xm';

%tm=tm';

 x=xn;

 t=tn;

%numericos

fn=double(subs(ffs,{xx,tt},{xn,tn}));

un=double(subs(us,{xx,tt},{xn,tn}));

 pn = double(subs(ps,{xx,tt},{xn,tn}));

%f=fn;

%plot(pn)

A = zeros(J*N,N*J);

%B=zeros(J*N,N*J);

%A = zeros(J+N,J+N);

% pasar a avector las matrices

uu =un(:);

fvector = fn(:);

%avector =pn(:);

 avector =double(pn);

h=1/(J-1);

 ateorica = pn;

 uteorica = double(un);

%% Matriz

% for i = 1:J*N-2

% %for i = 1:J+N

% if i == 1 % t < h/2 (0,h/2)

% A(i,i:i+1) = [((uu(i+1) - uu(i+2)))                     (5*uu(i+1) - 2*uu(i+2) - 3*uu(i))]/(h^2);

% elseif i == J+N % t > 1-h/2 (1-h/2,1)

% A(i,i-1:i) = [(5*uu(i-1) - 2*uu(i-2) - 3*uu(i))       (uu(i-1)-uu(i-2))]/(h^2);

% else % h/2 < t < 1-h/2 (h/2,1-h/2)

% A(i,i-1:i+1) = [(-uu(i) + uu(i-1))            (uu(i+1)-2*uu(i)+un(i-1))       (uu(i+1)-uu(i))]/(2*h^2);

% end

% end

%whos uu

%v=diag((double(-uu(2:J*N-2) ) +double(uu(1:(J*N-2)-1)))/(2*h^2),-1)

v1=diag(double([(-uu(2:J*N-2)+uu(1:(J*N-2)-1))/(2*h^2); 0;0]),-1);

v2=diag(double([0; (uu(3:(J*N-2)+1)-2*uu(2:(J*N-2)) + uu(1:(J*N-2)-1)  )/(2*h^2);0;0]));

v3=diag(double([0; (uu(3:(J*N-2)+1)-uu(2:(J*N-2)))/(2*h^2); 0]),1);

A=v1+v2+v3;

A(1,1:2)= [((uu(1+1) - uu(1+2)))                     (5*uu(1+1) - 2*uu(1+2) - 3*uu(1))]/(h^2);

A(J+N,J+N-1:J+N) = [(5*uu(J+N-1) - 2*uu(J+N-2) - 3*uu(J+N))       (uu(J+N-1)-uu(J+N-2))]/(h^2);

%B(J*N/2,J*N/2+1)=0;

% A

% m=A-B

% [f,c]=find(m~=0)

%V=[v; zeros(2,size(v,1))]

%whos m

aa = double(pinv(A)*fvector);

whos

%delete(gcp('nocreate'))

plot(aa)

function  [A,ateorica,uteorica,fvector,avector,x,t]=calorparabolica

close all

J=30;

%funcion a la perfeccion hasta J=30, pero puedo intentar con otro

%numeros mayores

N=J;

xn=linspace(0,1,N*J);

tn=linspace(0,1,J*N);

syms xx tt

us = vpa(exp(tt).*sin(pi.*xx),3);

ps = vpa(40*xx.^2./(40*xx.^2+1)+1+0.*tt,3);

ffs=vpa(-diff(ps.*diff(us,xx),xx)+diff(us,tt),3);

 x=xn;

 t=tn;

fn=double(subs(ffs,{xx,tt},{xn,tn}));

un=double(subs(us,{xx,tt},{xn,tn}));

 pn = double(subs(ps,{xx,tt},{xn,tn}));

A = zeros(J*N,N*J);

uu =un(:);

fvector = fn(:);

 avector =double(pn);

h=1/(J-1);

 ateorica = pn;

 uteorica = double(un);

v1=diag(double([(-uu(2:J*N-2)+uu(1:(J*N-2)-1))/(2*h^2); 0;0]),-1);

v2=diag(double([0; (uu(3:(J*N-2)+1)-2*uu(2:(J*N-2)) + uu(1:(J*N-2)-1)  )/(2*h^2);0;0]));

v3=diag(double([0; (uu(3:(J*N-2)+1)-uu(2:(J*N-2)))/(2*h^2); 0]),1);

A=v1+v2+v3;

A(1,1:2)= [((uu(1+1) - uu(1+2)))                     (5*uu(1+1) - 2*uu(1+2) - 3*uu(1))]/(h^2);

A(J+N,J+N-1:J+N) = [(5*uu(J+N-1) - 2*uu(J+N-2) - 3*uu(J+N))       (uu(J+N-1)-uu(J+N-2))]/(h^2);

aa = double(pinv(A)*fvector);

plot(aa)