////////////////////////////////////////////////////////////////////////////////
// Aproximação numérica para a solução estacionária da equação                //
// do calor em uma dimensão. A partir de operadores de diferença              //
// finita em n pontos uniformemente espaçados.                                //
//                                                                            //
// Vamos considerar a equação                                                 //
//                           (k(x)u'(x))'=-f(x),                              //
// k(x) é o coeficiente de difusividade  térmica e f é uma fonte              //
// (ou sorvedouro) de calor.                                                  //
//                                                                            //
// Aplicando a operação de derivação no lado direito da equação, nota-se que  // 
// ela pode ser reescrita como                                                //
//                      k'(x)u'(x)+k(x)u''(x)=-f(x).                          //
//                                                                            //
// O domínio da equação é o interior de uma barra com extremidades em 0 e L.  //
// Nas extremidades serão considerardas condições de contorno.                //
////////////////////////////////////////////////////////////////////////////////


// Comprimento da barra.
L=1;

// Número de pontos na aproximação. 
// A solução é um vetor de n componentes. u(1) corresponde à solução em x=0 e 
// u(n) corresponde à solução em x=L.
n=1000;


// Espaçamento h.
h=L/(n-1);


// Definição do coeficiente de difusividade térmica.
function z=k(x)
    z= 2-sin((%pi/L)*x);
endfunction

// Definição da(s) fonte(s) térmica(s).
function z=f(x)
    z=150*exp(-10*(x-0.6)^2)-0*exp(-10*(x-0.5)^2);
endfunction

// Operador de diferença finita para d(.)/dx:
// f'(x)=( f(x+h)-f(x-h) )/(2h) + O(h^2)

// Operador de diferença finita para d^2(.)/(dx^2):
// f''(x)=( f(x+h)-2f(x)+f(x-h) )/(h^2) + O(h^2)

// O operador é utilizado em pontos internos, ou seja, para i=2,3,...,n-1.
// Em i=1 e i=n utilizamos as condições de contorno.

// A matriz dos coeficientes
U=zeros(n,n);

// Condição de contorno -> u(0)=0 e u'(L)=0.  
U(1,1)=1; // em x=0.

U(n,n)=3;U(n,n-1)=-4;U(n,n-2)=1; // em x=L.

// Operador de diferença finita. 
for i=2:n-1
    U(i,i-1)=-(k((i+1)*h)-4*k(i*h)-k((i-1)*h));
    U(i,i)=-8*k(i*h);
    U(i,i+1)=(k((i+1)*h)+4*k(i*h)-k((i-1)*h));
end;

// Coluna das constantes

// Condição de contorno -> u(0)=0 e u'(L)=0   
b(1)=0; // em x=0 .
b(n)=0; // em x=L.

// Lado direito da equação (fonte térmica). 
for i=2:n-1
  b(i)=-(4*h^2)*f((i-1)*h)/k((i-1)*h);
end


// Solução numérica

u=U\b;


plot((0:n-1)*h,u);
xgrid

// 