//////////////////////////////////////////////////////////////////
//////// problemes aux limites pour l'équation de Poisson 1D /////////
//////////////////////////////////////////////////////////////////
//
//
//  PROBLEME 1:  CONDITIONS LIMITES DE DIRICHLET (méthode de tir et différences finies)
//
//
deff('Y=fode(t,X)','Y=[X(2);c(t)*X(1)-f(t)]');
//
deff('y=f(x)','y=(%pi^2+1)*sin(%pi*x)')
deff('y=c(x)','y=1*ones(x)');
deff('y=uex(x)','y=sin(%pi*x)')
//
// solution exacte
//
xset('window',1)
clf()
t=0:0.01:1;
clf()
plot2d(t,uex(t),-1)
//
//  METHODE 1: TIR
function u1=phi(lambda)
//
ulambda=ode([0;lambda],0,t,fode)
u1=ulambda(1,$)
endfunction
//
//
lambda0=fsolve(0,phi);
uapprox=ode([0;lambda0],0,t,fode)
//
plot2d(t,uapprox(1,:),2)
//
//  METHODE 2: DIFFERENCES FINIES
N=100;
h=1/(N+1);xi=h:h:N*h;
//
A=1/(h^2)*toeplitz([2,-1,zeros(1,N-2)])+diag(c(xi));
F=f(xi');
//
U=inv(A)*F;
plot2d([0;xi';1],[0;U;0],3)
//
legends (['exact','tir','differences finies CL Dirichlet'],[-1,2,3]);
//
//  PROBLEME 2:  CONDITIONS LIMITES MIXTES DIRICHLET NEUMAN (différences finies)
//
deff('y=c2(x)','y=x');
deff('y=f2(x)','y=(x+%pi^2/4*ones(x)).*sin(%pi/2*x)');
deff('y=uex2(x)','y=sin(%pi/2*x)');
N=100;
h=1/(N+1);xi=h:h:(N+1)*h;
//
A=1/(h^2)*toeplitz([2,-1,zeros(1,N-1)])+diag(c2(xi));A(N+1,N)=-1/h^2+A(N+1,N);
F=f2(xi');
//
xset('window',2)
clf()
t=0:0.01:1;
plot2d(t,uex2(t),-1)
U=inv(A)*F;
plot2d([0;xi'],[0;U],1)
//
legends (['differences finies CL mixtes','uexact'],[1,-1]);
//
//
//  PROBLEME 3:  CONDITIONS LIMITES NEUMAN (différences finies)
//
deff('y=c3(x)','y=1*ones(x)');
deff('y=f3(x)','y=-2*ones(x)+12*x-12*(x.^2)+(x.^2).*(1-x).^2');
deff('y=uex3(x)','y=(x.^2).*(1-x).^2');
N=100;
h=1/(N+1);xi=0:h:(N+1)*h;
//
A=1/(h^2)*toeplitz([2,-1,zeros(1,N)])+diag(c3(xi));
A(N+2,N+1)=-1/h^2+A(N+1,N);
A(1,2)=-1/h^2+A(1,2);
F=f3(xi');
//
xset('window',3)
clf()
t=0:0.01:1;
plot2d(t,uex3(t),-1)
U=inv(A)*F;
plot2d([xi'],[U],2)
//
legends (['differences finies CL neuman','exacte'],[2,-1]);
//











