Newton r

function[x,it]=Newton_R(x0,tol)
x=sym('x');
g=@(x)f1(x);
h(x)=diff(g(x),x,1);
xk=x0-f1(x0)/(vpa(h(x0),3));
i=1;
while (abs(xk-x0)>tol)
x0=xk;
xk=x0-f1(x0)/(vpa(h(x0),3));
i=i+1;
end
x=xk; it=i;
return

jacobi
function[X]=jacobi(A,b,kmax)
n=length(b); x0=zeros(n,1);
D=diag(diag(A)); E=-tril(A-D); F=-triu(A-D);
Ej=D\(E+F); qj=D\b;
i=1;
xk=Ej*x0+qj;
while i<=kmax
x0=xk;
xk=Ej*x0+qj;
i=i+1;
end
X=xk;
return

gauss s
function[X]=Gauss_Seidel(A,b,kmax)
n=length(b); x0=zeros(n,1);
D=diag(diag(A)); E=-tril(A-D); F=-triu(A-D);
EGS=(D-E)\F; qGS=(D-E)\b;
i=1;
xk=EGS*x0+qGS;
while i<=kmax
x0=xk;
xk=EGS*x0+qGS;
i=i+1;
end
X=xk;
return

gauss
function[U,f]=Gauss(A,b)
n=length(b);
for k=1:n-1
for i=(k+1):n
m=A(i,k)/A(k,k);
for j=k:n
A(i,j)=A(i,j)-m*A(k,j);
end
b(i)=b(i)-m*b(k);
end
end
U=A; f=b;
return

back substitution
function[X]=Bsubstitution(A,b)
n=length(b);
[U,f]=Gauss(A,b);
for i=n:-1:1
x(i)=f(i);
for j=(i+1):n
x(i)=x(i)-U(i,j)*x(j);
end
x(i)=x(i)/U(i,i);
end
X=x;
Return

Downsubt
function[X]=downsubstitution(c,b)
n=length(b);
%[U,f]=Gauss(A,b);
U=c; f=b;
for i=1:n
x(i)=f(i);
for j=1:(i-1)
x(i)=x(i)-U(i,j)*x(j);
end
x(i)=x(i)/U(i,i);
end
X=x;
return
polynome
function[L]=POL(x,z)
%x est le pont d interpolation
%z is a target points
n=length(x); m=length(z); L=[];
for k=1:m
T=[];
for p=1:n
v=x; a=v(p); v(p)=[];
A=[]; B=[];
for j=1:(n-1)
A=[A,z(k)-v(j)]; B=[B,a-v(j)];
end
l(p)=prod(A)/prod(B); T=[T;l(p)];
end
L=[L,T];
end
return

lagrange
function[P]=Lagrange(x,y,z)
%P: the polynomial evalueted at all targets z
L=POL(x,z);
P=L'*y;
%plot(x,P,"-r",x,y,"-b")
return

DD
function[b]=DD(x,y)
%cette fonction implemente l algo de dfference de division
%x:vecteur des point d interpolation
%y:valeur de vecteurs fonction
%b: solution des coefficient de Newton
n=length(x); b=y;
for i=2:n
for j=2:i
b(i)=(b(i)-b(j-1))/(x(i)-x(j-1));
end
end
return
newton
function[p]=Newton(x,y,z)
%this function implements the algoritm of divided
difference
%z:vector of target point
%p: the polynomial at all targets
b=DD(x,y);
n=length(b);
for i=1:length(z)
p(i)=b(n);
for k=(n-1):-1:1
p(i)=p(i)*(z(i)-x(k))+b(k);
end
end
return

F
function[x]=F(b,t,l,u)
S=0;
sk=0;
k=0;
for i=0:(b-2)
sk=sk+1;
k=i+1;
for k=1:t
for p=l:u
if S>0
S=S+(b^p)*(sk*(b^-k));
else
S=S-(b^p)*(sk*(b^-k));
end
end
end
end
x=S;
return

relaxation
function[X]=Relaxation(A,b,kmax,w)
n=length(b); x0=zeros(n,1);
D=diag(diag(A)); E=-tril(A-D); F=-triu(A-D);
ER=((1/w)*D-E)\((1/w-1)*D+F); qR=((1/w)*D-E)\b;
k=1;
xk=ER*x0+qR;
while k<=kmax
x0=xk;
xk=ER*x0+qR;
k=k+1;
end
X=xk;
return

bissection
function[x,it]=bissection(a,b,tol)
c=(a+b)/2; fc=f2(c); i=1;
while abs(b-a)>tol
if f2(a)*f2(c)<0
b=c;
c=(a+b)/2;
i=i+1;
else
a=c;
c=(a+b)/2;
i=i+1;
end
end
x=c; it=i;
return

secantes
function[x,it]=secante(x0,x1,tol)
xk=x1-(fe(x1))*(x1-x0)/(fe(x1)-fe(x0));
k=1;
while (abs(xk-x1)>tol)
x0=x1;
x1=xk;
xk=x1-(fe(x1))*(x1-x0)/(fe(x1)-fe(x0));
k=k+1;
end
x=xk; it=k;
return
point fixe
function[X,it]=point_fixe(g,x0,tol)
n=0;
x=x0;
xk=g(x);
k=1;
while abs(xk-x0)>tol
xk=x0;
xk=g(xk);
k=k+1;
n=n+1;
end
X=xk; it=n;
return

discriminent
function[x]=discriminant(a,b)
delta=a^2-4*b;
if delta<0
x='pas de solution';
else
x=[(-a+sqrt(delta))/2,(-a-sqrt(delta))/2];
end
return

suites
function[X]=suite(u0,kmax)
k=1; uk=sqrt(5*u0+2);
while (k<=kmax)
u0=uk;
uk=sqrt(5*u0+2);
k=k+1;
end
X=uk;
return

mat
function[A]=mat(a,b,c,N)
A=zeros(N,N);
for i=2:N
 A(i,i)=a;
A(i-1,i)=b;
A(i,i-1)=c;
end
A(1,1)=a;
return

f1
function[y]=f1(x)
s=6000;
for k=1:5
s=s-1000*(1+x).^k;
end
y=s;
return

fac
function[n]=fact(i)
n=1;
j=1;
while j<=i
n=n*j;
j=j+1;
end

hilbert
%exercise1
%question 1
%generons la matrice hilbert d ordre 10
H=hilb(10);
%vecteur diagonal
d=diag(H);
%decomposition traiangulaire [U,L]
L=tril(H); U=triu(H);
%affichage du resultat
disp("diagonal vector(d):");
disp(d);
disp("lower triangular matrix(L):");
disp(L);
disp("upper triangular matrix(U):");
disp(U);
%question 2
x=ones(10,1);
kmax=10000;
n=length(x);
b=H*x;
%affichage du resultat
disp("solution of(b):");
disp(b);
%auestion 3
%resolution avec pertubation sb
sb=[zeros(9,1);1e-4];
c=b+sb;
y=jacobi(H,c,kmax);
%affichage de resultat
disp("pertubetor solution(y):");
disp(y);

%donnons l estimation de la condition K(H)