Equa_diff1.m

function Equation_diff1(n)

x=[0:pi/n:pi]';             # vecteur contenant les coordonnées des differents noeuds

K=zeros(n+1); ;         # initialisation de la matrice de raideur K: n+1x n+1

F=zeros(n+1,1); ;       # initialisation de F: vecteur a n+1 lignes

 

for i=1:n ;                 # initialisation du compteur de boucle

j=i+1;

t=[i,j];             # matrice de connectivite

x1=x(i);          # coordonnées i de l'élément

x2=x(j);          # coordonnées j de l'élément

[Ke,Fe]=Calcul_Mat_element(x1,x2);         # calcul de Ke et de Fe

K(t,t)=K(t,t)+Ke;             # assemblage de la matrice de raideur

F(t)=F(t)+Fe;                 # assemblage du vecteur de charge

endfor; # fin de boucle

##################################### Application des conditions aux limites

F=F-1*K(:,1);           # prise en compte de la valeur u(0)=1

K(n+1,:)=[];             # suppression de la n+1 ieme ligne de K

K(:,n+1)=[];             # suppression de la n+1 ieme colonne de K

K(1,:)=[];                 # suppression de ligne 1 de K

K(:,1)=[];                 # suppression de colonne 1 de K

F(n+1)=[];               # suppression de la n+1 ieme ligne de F

F(1)=[];                   # suppression de la ligne 1 de K

#####################################

U=K\F;                    # calcul de U

U=[1;U;0];                # on complete U avec les valeurs aux bornes

plot(x,U);                 # tracée de U en fonction de x

endfunction

Calcul_Mat_element.m

function [Ke, Fe]=Calcul_Mat_element(x1, x2)

delta=x2-x1;

Ke1=1/(delta)*[1 , -1;-1, 1];

Ke2=1/2*[x2+2*x1, -x2-2*x1; 2*x2+x1, -2*x2-x1];

Ke3=-7/6*delta*[2,1;1,2];

Ke=Ke1+Ke2+Ke3;

A=-0.5*delta*sin(2*x1)+0.5*sin(delta)*sin(x2+x1);

B=0.5*delta*sin(2*x2)-0.5*sin(delta)*sin(x2+x1);

Fe=-1/(delta)*[A; B];

return

endfunction