Esercizi svolti in MATLAB, Calcolo numerico

C=L';

for i=n:-1:1

x(i)=y(i);

for j =i+1:n

x(i)=x(i)-C(i,j)*x(j);

end

x(i)=x(i)/C(i,i);

end

x=x';

end

choleky_test

(permette di richiamare la function dal command window.Basta scrivere il nome

della function "[ L,x ] = CHOLO( A,b,n )" e premere INVIO. E' anche un modo di

verifica della function)

A=[ 4 2 -2;2 10 -7;-2 -7 9]

b=[-8 -7 2]';

[ L,x ] = CHOLO( A,b,n )

--------------------------------------------------------------

function [ C,p ] = DIFF_DIVISE( x,y,xx,n )

%DIFF_DIVISE Summary of this function goes here

% Detailed explanation goes here

C=zeros(n,n);

for i=1:n

C(i,1)=y(i);

end

for j=2:n

for i=j:n

C(i,j)=(C(i,j-1)-C(i-1,j-1))/(x(i)-x(i-j+1));

end

end

m=length(xx);

p=ones(m,1)*C(n,n)

for j=n-1:-1:1

p=p.*(xx-x(j))+C(j,j);

end

end

DIFF_DIVISE _TEST

x=[-1 2 -2 3 ]';

y=[-1 2 -2 4 ]';

xx=[0 1 2 4 5 6 0.1 ]';

n=4;

[ C,p ] = DIFF_DIVISE( x,y,xx,n )

-----------------------------------------------------------------

function [ C,yy] = diffdivise( x,y,xx )

%DIFFDIVISE Summary of this function goes here

% Detailed explanation goes here

n=length(x);

C=zeros(n,1);

C(:,1)=y(:);

for j=2:n

for i=j:n

C(i,j)=(C(i,j-1)- C(i-1,j-1))/(x(i)-x(i-j+1));

end

end

%HORNER

m=length(xx);

yy=zeros(1,m);

for i=1:m

p=C(1,1);

for j=2:n

p=p+C(j,j)*prod(xx-x(1:j-1));

end

yy(i)=p;

end

end

test

x=[0 1 2 3 ]

y=[5 -2 4 18]

xx=0.5

[ C,yy] = diffdivise( x,y,xx )

----------------------------------------------------------------

function [ x,y ] = EULERO(f,x0,y0,b,h )

%EULERO Summary of this function goes here

% Detailed explanation goes here

m=(b-x0)/h;

x=0;

y(1)=y0;

for n=2:m+1

x(n)=x0+(n-1)*h;

y(n)=y(n-1)+h*feval(f,x(n-1),y(n-1));

end

x=x';

y=y';

T=[x y]

end

test

f=inline('-y*exp(x)');

b=0.2;

x0=0;

y0=1

h=0.1;

[ x,y ] = EULERO(f,x0,y0,b,h )

----------------------------------------------------------------

function [x,y ] = EURLERO_MODIFIC(f,x0,y0,b,h )

%EURLERO_MODIFIC Summary of this function goes here

% Detailed explanation goes here

x=x0;

y(1)=0;

m=(b-x0)/h;

for n=2:m+1

k1=feval(f,x(n-1),y(n-1));

k2=feval(f,x(n-1) + h*0.5,y(n-1)+h*0.5*k1);

y(n)=y(n-1)+h*k2;

x(n)=x0+(n-1)*h;

end

x=x';

y=y';

T=[x, y]

end

test

f=inline('4*x.^3*sqrt(1-y.^2)');

x0=0;

y0=0;

b=0.2;

h=0.1;

[x,y ] = EURLERO_MODIFIC(f,x0,y0,b,h )

----------------------------------------------------------------

function [L Uy ] = GAUS_PIVOTING(Ab,n )

%GAUS_PIVOTING Summary of this function goes here

% Detailed explanation goes here

n=3;

L=ones(n,n);

m=length(Ab);

Uy=Ab;

% ALGORITMO DI GAUSS CON PIVOTING

for k=1:n-1

[pivot indice]=max(abs(Ab(k:n,k)));

riga=indice+k-1;

if riga~= k

Ab([riga k],:)=Ab([k riga],:);

end

for i=k+1:n

Ab(i,k)=Ab(i,k)/Ab(k,k);

for j=k+1:n+1

Ab(i,j)=Ab(i,j)-Ab(i,k)*Ab(k,j);

end

Uy(i,j)=Ab(i,j);

end

end

% COSTRUZIONE L

for i=1:n

L(i,i)=1;

for j=i+1:n

L(i,j)=0;

L(j,i)=Ab(j,i);

Ab(j,i)=0;

end

end

% COSTRUZIONE Uy e U

for k=1:n-1

for i=k+1:n

R(i,k)=Ab(i,k)/Ab(k,k);

for j=k:n+1

Ab(i,j)=Ab(i,j)-R(i,k)*Ab(k,j);

Uy(i,j)=Ab(i,j);

end

end

end

U=Uy(1:n, 1:m-1)

end

test

A=[5 5 3;2 2 -2 ;3 6 5 ];

b=[1 1 1 ]';

Ab=[A b];

n=3;

[L Uy ] = GAUS_PIVOTING(Ab,n )

-------------------------------------------------------

function [A ] = gauss_pivoting( A,b )

%GAUSS_PIVOTING Summary of this function goes here

% Detailed explanation goes here

n=length(A);

x=zeros(n,1);

[m, n]=size(A);

for k=1:n

[pivot indice]=max(abs(A(k:n,k)));

riga=indice+k-1;

if riga~= k

A([riga k],:)=A([k riga],:);

b([riga k])=b([k riga]);

end

end

for i=k+1:n

A(i,k)=A(i,k)/A(k,k);

for j=k+1:n

A(i,j)=A(i,j)-A(k,j)*A(i,k);

end

b(i)=b(i)-b(k)*A(i,k)

end

end

test

A=[1 -1 2;3 0 1;2 2 2]

b=[1 1 1 ]

[ A ] = gauss_pivoting( A,b )

-------------------------------------------------------

function [ L,U,x ] = gauss( A,b )

%FATTORIZZO MATRICE A=L*U

% E L'EQUIVALENTE DI GAUSS SENZA PIVOTING

n=length(A);

L=zeros(n,n);

U=zeros(n,n);

for j=1:n

for i=1:j

U(i,j)= A(i,j)-L(i,1:i-1)*U(1:i-1,j);

L(i,i)=1;

for i=j+1:n

L(i,j)=(A(i,j)-L(i,1:j-1)*U(1:j-1,j))/U(j,j);

end

end

end

% RISOLVO IL SISTEMA (Ly=b Ux=y)

y=zeros(n,1);

x=zeros(n,1);

for i=1:n

y(i)=b(1)/L(i,i);

end

for i=2:n k=1:i-1;

y(i)=(b(i)-L(i,k)*y(k))/L(i,i);

end

for i=n:-1:1

x(i)=y(i)/U(i,i);

end

for i=n-1:-1:1 k=i+1:n;

x(i)=(y(i)-U(i,k)*x(k))/U(i,i);

end

Uy=[U y]

end

test

A=[ 1 0 2;2 2 6;3 -4 1];

b=[0 0 1 ];

[ L,U,x] = gauss( A,b )

-------------------------------------------------------

function [x,r,k ] = GRAD_CONIUG(A,b,x0,alfa,toll,kmax )

%GRAD_CONIUG Summary of this function goes here

% Detailed explanation goes here

n=length(A);

k=0;

xk=x0;

test=norm(b-A*xk);

toll=toll*norm(b);

rk=b-A*xk

x=[];

r=[];

while test > toll && k<kmax

k=k+1;

xk=xk+alfa*rk;

rk=b-A*xk;

x=[x;xk];

r=[r;rk];

end

end

test

n=30;

v=ones(n,1);

A=spdiags ([-0.2*v 4*v -0.2*v ], -1:1,n,n) ;

b=ones(n,1);

toll=0.0000001

kmax=3;

x0=zeros(n,1);

alfa=1.3;

[x,r,k ] = GRAD_CONIUG(A,b,x0,alfa,toll,kmax )

--------------------------------------------------------

function [ x,y ] = HEUN( f,x0,y0,b,h )

%HEUN Summary of this function goes here

% Detailed explanation goes here

m=(b-x0)/h;

x0=0;

y(1)=y0;

for n=2:m+1

x(n)=x0+(n-1)*h;

k1=feval(f,x(n-1),y(n-1));

k2=feval(f,x(n),y(n-1)+h*k1);

y(n)=y(n-1)+h*(k1+k2)*0.5;

end

end

test

f=inline('4*x.^3*sqrt(1-y.^2) ');

x0=0;

b=0.2;

y0=0;

h=0.1;

[ x,y ] = HEUN( f,x0,y0,b,h )

------------------------------------------------------------

function [ x,k ] = JACOB(A,b,x0,n,toll,kmax )

%JACOBI_ Summary of this function goes here

% Detailed explanation goes here

n=length(A);

k=0;

x=x0;

test=toll+1;

while test>toll && k< kmax

k=k+1;

for i=1:n

xstar(i)=0;

for j=1:i-1

xstar(i)=xstar(i)+A(i,j)*x(j);

end

for j=i+1:n

xstar(i)=xstar(i)+A(i,j)*x(j);

end

xstar(i)=(b(i)-xstar(i))/A(i,i);

end

test=norm(xstar(i)-x);

x=xstar';

end

end

------------------------------------------------------

function [ yy ] = lagrange(x,y,xx)

%LAGRANGE valuta il polinomio in un insieme di punti xx

n=length(x);

m=length(xx);

yy = zeros(m,1);

%calcolare il polinomio generico ci Lagrange

for i=1:n

for j=1:n

if j ~= i

yy =yy+y(i)*(xx-x(j))/(x(i)-x(j));

end

end

end

----------------------------------------------------------------

function [ T,yy ] = NEVILLE(x,y,xx,n )

%NEVILLE Summary of this function goes here

% Detailed explanation goes here

n=length(x);

T=zeros(n,n);

for i=1:n

T(i,1)=y(i);

end

for k=1:n

for i=1:n-k

T(i,k+1)=T(i+1,k)-(T(i+1,k)-T(i,k))*(x(i+k)-xx)/(x(i+k)-x(i));

end

end

yy=T(1,n);

end

test

x=[3 1 2]';

y=[10 2 1 ]';

n=3;

xx= 3/2;

[ T,yy ] = NEVILLE(x,y,xx,n )

----------------------------------------------------------------

function [ p ] = newLAGRANGE( x,y,a)

%NEWLAGRANGE Summary of this function goes here

% Detailed explanation goes here

n=length(x);

m=length(a)

p(1,:)=zeros(1,m);

for k=1:m

for i=1:n

L(i)=1

for j=1:i-1

L(i)=L(i)*(a(k)-x(j))/(x(i)-x(j));

end

for j=i+1:n

L(i)=L(i)*(a(k)-x(j))/(x(i)-x(j));

end

test

x=[-2 0 1 2 ]

y=[4 10 10 16]

a=[1 2 3 4 5 ]

[ p ] = newLAGRANGE( x,y,a)

-----------------------------------------------------------

function [ x,n,res ] = NEWTO( f,df,x0,toll,nmax )

%NEWTO Summary of this function goes here

% Detailed explanation goes here

n=0;

diff=toll+1

x(1)=x0;

while diff > toll && n < nmax

n=n+1;

fx=feval(f,x(n));

dfx=feval(df,x(n));

diff=-fx/dfx;

x(n+1)=x(n)+diff;

res=feval(f,x(n+1));

x=x';

diff=abs(diff);

end

end

test

f=inline('sin(x)+x-pi');

df=inline('cos(x)+1');

x0=2;

toll=0.0000001;

nmax=3;

[ x,n,res ] = NEWTO( f,df,x0,toll,nmax )

----------------------------------------------------------

function [ xv resv ] = newton(f,df,x0,toll,nmax )

x0=input('dammi il valore iniziale x0= ');

toll=input('dammi il valore tolleranza toll= ');

nmax=input('dammi il valore iterate nmax= ');

f=inline('x^3-2*x+2');

df=inline('3*x^2-2');

diff=toll+1;

xv=[];

resv=[];

n=-1;

x=x0;

while diff>=toll && n < nmax

n = n+1; % incrementa il contatore delle iterate

fx=feval(f,x);

dfx=feval(df,x);

x = x -(fx/dfx); % Calcola il valore x_{n+1}

xv=[xv x];

fx = feval(f,x); % Calcola il valore f(x_n)

resv =[resv fx];

diff=abs(fx/dfx);

----------------------------------------------------------

function [ xv resv ] = puntofisso(f,g,x0,toll,nmax )

%PUNTOFISSO Summary of this function goes here

%

xv=[];

resv=[];

res=toll+1;

n=-1;

x=x0;

while res>toll && n< nmax

n=n+1;

xold=x;

gx=feval(g,xold);

x=gx;

xv=[xv;x];

fxnew=feval(f,x);

resv=[resv;fxnew];

res=abs(xold-x);

end

test

f=inline('exp(-3*x)-x+1');

g=inline('(exp(-3*x)*(3*x+1)+1)/(3*exp(-3*x)+1)');

x0=1;

toll=0.0000000000000005;

nmax=5;

[ xv resv ] = puntofisso(f,g,x0,toll,nmax )

-----------------------------------------------------------

function [ a ] = REGRESSION(x,y,w )

%REGRESSION Summary of this function goes here

% Detailed explanation goes here

n=length(x);

s=sqrt(w);

A=[s s.*x s.*x.*x]

b=[s.*y];

C=A'*A;

D=A'*b;

%RISOLVO IL SISTEMA Ca=D;

a=inv(C)*D;

a=a';

end

test

x=[ 0 1/2 1 3/2 2 ]'

y=[0 1/2 1 0 -1 ]'

w=[1 16 1 16 1 ]'

s=sqrt(w)

[ a ] = REGRESSION(x,y,w)

plot(x,y,'r-')

------------------------------------------------------------

function [ Istar T ] = ROMBER(f,a,b,n,m )

%ROMBER Summary of this function