%% Práctica 2. Ecuaciones en Derivadas Parciales.
% Veremos cómo implementar métodos numéricos basados en diferencias
% finitas para resolver aproximadamente la ecuación de ondas, 
% la ecuación del calor y la ecuación de Poisson.

%% Método explícito para el problema mixto para la ecuación de ondas
% En esta sección vamos a implementar en Matlab/Octave el método explicado en 
% teoría. Debemos tener en cuenta que aunque matemáticamente se pueda definir un 
% subíndice de una variable como cero, los elementos de las matrices (y los índices 
% en Matlab/Octave) empiezan en 1.

%%
% Como ejemplo, nos planteamos el problema
%
% (parcial segunda de u respecto de t dos veces) - 4*(parcial segunda de u
% respecto de x dos veces) = 0 para t>=0, x en [0,6],
%
% u(0,x)=x*(x-6)^2,
%
% parcial de u respecto de t en (0,x)=0, para x en [0,6],
%
% u(t,0)=u(t,2)=0, para t>=0.
% donde discretizamos en el dominio [0,2] x [0,6] de forma que en el dominio del
% tiempo haya 2 subintervalos y en el dominio de x haya 3 subintervalos.

% Comenzamos definiendo la función phi que define la deformación inicial de la cuerda.
% La definimos de forma que se pueda usar en vectores.
phi=@(x) x.*(x-6).^2;

%%
% La solución del problema será una función u(t,x) donde la variable t está en el
% intervalo [0,T] dividido en N intervalos y la variable x está en el
% intervalo [0,ele] dividido en M intervalos.
T=2;
ele=6;
N=2;
M=3;

%%
% Como el intervalo [0,T] está dividido en N subintervalos, cada uno de
% estos subintervalos tendrá longitud Deltat=T/N. De la misma manera, 
% como el intervalo [0,ele] está dividido en M subintervalos, cada uno tendrá
% longitud Deltax=ele/M. Con esto podemos definir los puntos que utilizamos
% para dividir los intervalos (y que dan la 'red uniforme' de la que se
% habla en la teoría).
Deltat=T/N;
Deltax=ele/M;
nodost=[0:Deltat:T];
nodosx=[0:Deltax:ele];

%%
% Definimos el coeficiente c a partir de la ecuación en derivadas parciales y el
% valor de  q que se usará en las ecuaciones de diferencias finitas.
c=2; % c^2 es el coeficiente de la parcial segunda de u respecto de x dos veces.
q=c*Deltat/Deltax;

%%
% Inicializamos una matriz en la que almacenaremos las aproximaciones de
% u en los puntos de la red uniforme.
U=zeros(N+1,M+1);

%% Condiciones iniciales para la ecuación de ondas.
% Definimos la deformación inicial de la cuerda. U(1,j)=phi( Deltax*(j-1) ) .
% Ahorramos un bucle aplicando phi a los nodos en el eje x.
U(1,:)=phi(nodosx);
% Ahora imponemos que la cuerda está sujeta en el extremo izquierdo x=0,
% que es equivalente a definir en un bucle U[i,1]=0 para i=1,2,...,N+1 . 
U(:,1)=0;
% También la cuerda está sujeta en el extremo derecho x=ele ,
% que es equivalente a definir en un bucle U[i,M+1]=0 para i=1,2,...,N+1 .
U(:,M+1)=0;

%% Primeros valores de U.
% Aquí usamos diferencias finitas para calcular U en el siguiente
% instante tras la deformación inicial de la cuerda.
% En la teoría estos valores se notaban con subíndices U1,j con i=1, 
% pero aquí usamos la notación U(2,j).
% No llegamos hasta U(2,M+1) porque ya sabemos que vale cero.
for j=2:M
    U(2,j)=(1-q^2)*U(1,j)+(q^2/2)*( U(1,j+1)+U(1,j-1) );
end

%% Bucle general para la ecuación de ondas.
% Ahora usamos el esquema de diferencias finitas visto en teoría 
% para el resto de puntos de la red uniforme.
% Empezamos calculando U(3,j), porque el U(2,j) ya se calculó en el
% anterior bucle.
for i=2:N
 for j=2:M  
 U(i+1,j)=2*(1-q^2)*U(i,j)+q^2*( U(i,j+1)+U(i,j-1) )-U(i-1,j) ;
 end
end

%% Resultados para la ecuación de ondas.
% La matriz U contiene los valores aproximados de la función u.
U
plot(nodosx, U(1,:) ) % Gráfica de la posición inicial i=1, t=0 .
disp('Pulse una tecla para ver la siguiente gráfica');
pause
plot(nodosx, U(2,:) )  %Gráfica del momento i=2, t=Deltax .
% Representar U.
% stem3(nodost,nodosx,U')

%%
% Ej. Se quiere hallar una aproximación a la solución del siguiente problema
% mixto para la ecuación de ondas
%
% (parcial segunda de u respecto de t dos veces) - (parcial segunda de u
% respecto de x dos veces) = 0 para t>=0, x en [0,2],
%
% u(0,x)=sen(pi*x/2),
%
% parcial de u respecto de t en (0,x)=0, para x en [0,2],
%
% u(t,0)=u(t,2)=0, para t>=0.
%
% Aplique el método explícito en diferencias finitas, dividiendo el
% dominio [0, 6] x [0, 2] en intervalos de longitud 2 en la variable t e
% intervalos de longitud 1 en el dominio de la variable x .
% (Sol. U=[ 0  1  0 ; 0  -3  0; 0  17  0; 0  -99  0 ] )

%% Método explícito para el problema mixto para la ecuación del calor
% Es este apartado nos planteamos hallar la solución aproximada de la
% ecuación del calor.
% Usaremos la misma notación y la misma discretización del dominio 
% que en la ecuación anterior. El cambio en la ecuación en derivadas parciales
% afectará a las ecuaciones en diferencias finitas.
%
% Nos planteamos el problema
% (parcial de u respecto a t) - (parcial segunda de u respecto a x dos
% veces) = 0, para t>=0, x en [0,3],
%
% u(0,x)=sen(pi*x/3), para x en [0,3],
%
% u(t,0)=u(t,3)=0, t>=0,
% donde consideramos el dominio del tiempo [0,9] dividido en 3
% subintervalos y el dominio [0,3] para la variable x, dividido en 3
% intervalos.

%% Inicialización para la ecuación del calor.
% Comenzamos definiendo la función phi que es la temperatura inicial de la
% barra, y discretizamos en los ejes del tiempo y del espacio.
phi=@(x) sin(pi*x/3);
T=9;
ele=3;
N=3;
M=3;
Deltat=T/N;
Deltax=ele/M;
nodost=[0:Deltat:T];
nodosx=[0:Deltax:ele];

%%
% Definimos el coeficiente de la ecuación en derivadas parciales y una
% matriz que contendrá las aproximaciones a los valores de u.
c=1; % Coeficiente de la parcial segunda de u respecto de x dos veces.
q=c*Deltat/(Deltax^2);
U=zeros(N+1,M+1);

%%
% Imponemos las condiciones de contorno.
U(1,:)=phi(nodosx);
U(:,1)=0;  %U(i,1)=0 para i=1,2,...,N+1 .
U(:,M+1)=0;  %U(i,M+1)=0 para i=1,2,...,N+1 . 

%%
% Ahora escribimos el bucle general usando el esquema en diferencias 
% finitas.
% En este caso tenemos que empezar por i=1, para calcular U(2,j). 
% Esto cambia con respecto al caso anterior.
for i=1:N
 for j=2:M  
 U(i+1,j)=(1-2*q)*U(i,j)+q*( U(i,j+1)+U(i,j-1) ) ;
 end
end

%% Resultados de la ecuación del calor.
% Para comprobar los resultados recordamos que
% -sqrt(3) = -1.7321, 
% 2*sqrt(3) = 3.4641, 
%  -4*sqrt(3)= -6.9282 .
U

%%
% Ej. Considere el siguiente problema mixto para la ecuación del calor
% (parcial de u respecto a t) - 2*(parcial segunda de u respecto a x dos
% veces) = 0, para t>=0, x en [0,4],
%
% u(0,x)=x*(4-x), para x en [0,4],
%
% u(t,0)=u(t,4)=0, t>=0.
%
% Halle una aproximación a la solución usando el método explícito en
% diferencias finitas en el dominio [0,4] x [0,4], dividiendo en intervalos
% de longitud 2 en el eje de variable t y en intervalos de longitud 1 en el
% eje de la variable x .
% (Sol. U=[ 0  3  4  3  0 ; 0  -5  -4  -5  0 ; 0  19  -12  19  0 ] )

%% Método para la ecuación de Poisson.
% En esta sección nos planteamos la ecuación de Poisson (laplaciano de u
% igual a una función dada f).
% Al implementar este método, tendremos que resolver un sistema
% de ecuaciones lineales para obtener las aproximaciones de la función U.
% Al plantear el sistema lineal, el vector término independiente contiene
% valores de f(xi,yj) y el vector de soluciones contiene valores de u.
% Será necesario reordenar los elementos de una matriz para que tengan
% forma de vector.

%%
% Nos planteamos el problema de aproximar una función u(x,y) definida en 
% un dominio D=[0,3]x[0,3] de forma que
%
% su laplaciano cumple
% (parcial segunda de u respecto a x dos veces) + (parcial segunda de u
% respecto a y dos veces) = 6(x+y), para (x,y) en D.
%
% u(x,y)= x^3+y^3, para (x,y) en la frontera de D.
% Para el método consideramos el dominio de la x dividido en tres
% subintervalos y el dominio de la y dividido en tres subintervalos.
  
%% Inicialización.
% Extremos del intervalo para el eje x y subintervalos.
a=0;
b=3;
N=3;
Deltax=(b-a)/N;
nodosx=[a:Deltax:b];

%%
% Extremos del intervalo para el eje y y subintervalos.
c=0;
d=3;
M=3;
Deltay=(d-c)/M;
nodosy=[c:Deltay:d];

% Definición de la función en la frontera.
phi=@(x,y)  x^3+y^3;

% Definición de la función f igual al laplaciano.
f=@(x,y) 6*(x+y);

% Inicializamos U
U=zeros(N+1,M+1);

% Inicializamos A y B para el sistema de ecuaciones.
A=zeros( (N+1)*(M+1), (N+1)*(M+1) );
B=zeros( (N+1)*(M+1), 1 );

%%
% La función q transforma los valores i=1,..,N+1, j=1,...,M+1 en valores de
% q=1,...,(N+1)*(M+1).
q=@(i,j) (j-1)*(N+1)+i;

%% Variables correspondientes a los puntos del borde.
% Usamos la función phi de las condiciones de contorno para definir
% las ecuaciones correspondientes a los valores de u en el borde.
for j=1:M+1
    % U(1,j)=phi(a, nodosy(j) );
    A( q(1,j),q(1,j) )=1;
    B( q(1,j) )=phi(a, nodosy(j) );
    
    % U(N+1,j)=phi(b, nodosy(j) );
    A( q(N+1,j),q(N+1,j) )=1;
    B( q(N+1,j) )=phi(b, nodosy(j) );
end

for i=1:N+1
    %U(i,1)=phi( nodosx(i), c);
    A( q(i,1),q(i,1) )=1;
    B( q(i,1) )=phi( nodosx(i), c);
    
    %U(i,M+1)=phi( nodosx(i), d);
    A( q(i,M+1),q(i,M+1) )=1;
    B( q(i,M+1) )=phi( nodosx(i), d);
end

%% Esquema de diferencias finitas.
% Montamos el sistema de ecuaciones que nos permite hallar los valores de
% U. Estas ecuaciones se deducen al aplicar diferencias finitas a los
% puntos interiores del dominio.

%%
% Queremos escribir la siguiente ecuación en forma matricial
%  (1/Deltax^2)*U(i+1,j) - (2/Deltax^2)*U(i,j) + (1/Deltax^2)*U(i-1,j) +
%  + (1/Deltay^2)*U(i,j+1)-(2/Deltay^2)*U(i,j)+(1/Deltay^2)*U(i,j-1) 
%  = f( nodosx( i ), nodosy( j ) )
% Cada valor de (i,j) nos da una ecuación distinta, que numeramos como
% ec=q(i,j) con término independiente  f( nodosx( i ), nodosy( j ) ) 
% y con coeficicientes (1/Deltax^2), ... que introducimos en la matriz A.

for i=2:N
    for j=2:M
    % Estamos en la ecuación con término independiente
    % f( nodosx( i ), nodosy( j ) )
    ec = q(i,j); %Número de la ecuación.      
    %
    A( ec, q(i+1,j) ) = (1/Deltax^2);
    A( ec, q(i,j) )      = -(2/Deltax^2) -(2/Deltay^2);
    A( ec, q(i-1,j) )  =  (1/Deltax^2);
    %    
    A( ec, q(i,j+1) ) = (1/Deltay^2);
    A( ec, q(i,j-1) )  = (1/Deltay^2);
    %    
    B( ec ) = f( nodosx( i ), nodosy( j ) );
    end
end

%%
% Resolvemos el sistema con la 'división izquierda'.
X=A\B

%%
% El vector X tiene (N+1)*(M+1) componentes que son 
% aproximaciones de los valores de U en los nodos.
% Recuperamos la matriz U.
for i=1:N+1
    for j=1:M+1
        U(i,j)=X( q(i,j) );
    end
end

U

%%
% Ej. Aplique el método en diferencias finitas para hallar una
% aproximación a la solución del siguiente problema de contorno 
% para la ecuación de Laplace en el dominio D formado por los puntos
% (x,y) que cumplen 0<x<4, 0<y<2 .
%
% (parcial segunda de u respecto a x dos veces) + (parcial segunda de u
% respecto a y dos veces) = 0, para (x,y) en D.
%
% u(x,y)=y^2-x^2, para (x,y) en la frontera de D.
%
% Para el método, divida el dominio en ambas variables en subintervalos
% de longitud 1.
% Compare los valores aproximados obtenidos con los de la función
% u(x,y)=y^2-x^2.
% (Sol. U = [ 0  1  4 ; -1     0     3 ; -4  -3  0 ; -9  -8  -5 ; -16  -15  -12
%  ] = valores_exactos ).

% Ecuaciones Diferenciales y Cálculo Numérico. (A. Palomares)
% Grado en Ingeniería de Tecnologías de Telecomunicación.
disp(date)