Esercizi svolti in MATLAB, Calcolo numerico

Codici MATLAB

Bisezione

function [ xn resn ] = bisezione(f,a,b,toll,nmax )

a=input('dammi il valore estremo sinistro');
b=input('dammi il valore estremo destro');
toll=input('dammi il valore tolleranza');
nmax=input('dammi il valore max iterazioni');
f=inline('x^3-x-5');
n=-1;
amp=toll+1;
xn=[];
resn=[];

while amp>=toll && n < nmax
n=n+1;
fa=feval(f,a);
fb=feval(f,b);
amp=abs(b-a);
x=a+(amp)*0.5;
fx=feval(f,x);
xn=[xn;x];
resn=[resn;fx];

if fa*fx>0
a=x ;
elseif fa*fx<0
b=x;
end

% Bisezione Summary of this function goes here
end
plot(xn,resn,'r*')

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

Bisez

function [ x res ] = bisez(f,a,b,toll,nmax )

a=input('dammi il valore di estremo sinitro a= ');
b=input('dammi il valore di estremo destro b= ');
toll=input('dammi il valore di toll= ');
nmax=input('dammi il valore di iterazioni nmax= ');
f=inline('x^3-x-5');
n=0;
chk=toll+1 ;
fa=feval(f,a);

while chk>=toll && n < nmax
n=n+1;
x(n)=a+(abs(b-a))*0.5 ;
fx=feval(f,x(n));
res(n)=fx;

if fa*fx>0
a=x(n);
elseif fa*fx<0
b=x(n);
end

% Bisez Summary of this function goes here
% Detailed explanation goes here

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

Cholo

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

% Cholo Summary of this function goes here
% Detailed explanation goes here
n=length(A);
L=zeros(n,n);

if A(1,1)<0
disp('no solution');
return
else
L(1,1)=sqrt(A(1,1));
end

for i=2:n
L(i,1)=A(i,1)/L(1,1);
end

for j=2:n
L(j,j)=A(j,j);
for k=1:j-1
L(j,j)=L(j,j)-L(j,k)*L(j,k);
end

if L(j,j)<0
disp('no solution');
return
else
L(j,j)=sqrt(L(j,j));
end

for i=j+1:n
L(i,j)=A(i,j);
for k=1:j-1
L(i,j)=L(i,j)-L(i,k)*L(j,k);
end
L(i,j)=L(i,j)/L(j,j);
end
end

% Risolvo il sistema Ly=b L'x=y

for i=1:n
y(i)=b(i);
for j=1:i-1
y(i)=y(i)-L(i,j)*y(j);
end
y(i)=y(i)/L(i,i);
end

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 )

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

Diff divise

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 )

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

Diffdivise

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 )

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

Eulero

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 )

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

Eurlero modific

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 )

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

Gaus p

function [L Uy ] = GAUS_P