Next: Risoluzione di sistemi
lineari e fattorizzazioni di matrici (pag.1) Up:Introduzione
Previous: Function
da 21 a 30
Intro Gen: Introduzione
Generale Home: Home
page
% METODO DI PIU' RIPIDA DISCESA PER LA RISOLUZIONE DI SISTEMI
LINEARI
function x=meripdis(A,b,n,e);
% x
vettore soluzione del sistema
% A vettore
dei coefficienti del sistema
% b
vettore dei termini noti
% n
dimensione di A
% e
é l'errore massimo che si deve commettere
for
i=1:n
somma=0;
for j=1:n
% determinazione di un intervallo
contenente
if j~=i
%
tutti gli autovalori della matrice
somma=somma+abs(A(i,j));
end
end
p1(i)=A(i,i)-somma;
p2(i)=A(i,i)+somma;
end
m=min(p1);
M=max(p2);
k=log(e*sqrt(m/M))/log((M-m)/(M+m));
x=diag(eye(n));
for
ki=1:k
r=b-A*x;
x=x+(((norm(r))^2)/(r'*A*r))*r;
end
x=x';
% METODO DEI TRAPEZI PER LA RISOLUZIONE DELLE EQUAZIONI
DIFFERENZIALI
function [y,risp]=mettrap(F,G,t,y,n,h,e);
% y
vettore soluzione dell'equazione differenziale
% risp
controlla se il metodo é applicabile
% F
é la funzione a due variabili F(t,y)
% G
é la derivata di F
% t
é la variabile indipendente della funzione y
% n
é il numero di suddivisioni dell'intervallo [a,b] fissato
% h
é l'ampiezza delle n suddivisioni
% e
é l'errore massimo che si deve commettere
risp=0;
wf1=feval(F,t(1),y(1));
i=2;
while i<=n+1 & risp==0
sol1=y(i-1);
wd=feval(G,t(i),sol1);
if
wd==2/h
risp=1;
else
wf2=feval(F,t(i),sol1);
wf3=sol1-y(i-1)-(h/2)*(wf1+wf2);
sol2=sol1-(wf3/(1-(h/2)*wd));
while risp==0 &
(abs(sol2-sol1))>=e | (abs(wf3)>=e)
sol1=sol2;
wf2=feval(F,t(i),sol1);
wf3=sol1-y(i-1)-(h/2)*(wf1+wf2);
sol2=sol1-(wf3/(1-(h/2)*wd));
end
wf1=wf2;
y(i)=sol2;
end
i=i+1;
end
% PROBLEMA LINEARE AI MINIMI QUADRATI
function [x,a,risp]=minquad(A,b,m,n);
% x
vettore soluzione del sistema
% a
norma 2 del residuo
% risp controlla
se il metodo é applicabile
% A
matrice dei coefficienti del sistema
% b
vettore dei termini noti
% m,n dimensioni
di A
[Q,R,risp]=fattqr(A,m,n);
% Fattorizzazione QR
if
risp==0
c=Q'*b;
x(n)=c(n)/R(n,n);
% Risoluzione del sistema R * x = c1 :
for
i=n-1:-1:1
a=0;
for j=i+1:n
a=a+R(i,j)*x(j);
end
x(i)=(c(i)-a)/R(i,i);
end
a=0;
for i=n+1:m
% Modulo del vettore residuo
a=a+(c(i))^2;
end
a=sqrt(a);
else
a=0;
x=0;
end
% SPLINE CUBICA COMPLETA INTERPOLANTE
function
[s,risp]=splcubco(x,y,xi,n,m1,k0,k1);
% s
vettore soluzione dei punti da interpolare
% risp controlla
se il metodo é applicabile
% x
vettore dei nodi interpolanti
% y
valori della funzione nei nodi x
% xi
vettore dei punti in cui calcolare la funzione interpolante
% n
dimensione di x , y
% m1 dimensione
di xi
% k0
valore della derivata prima nel nodo più piccolo
% k0
valore della derivata prima nel nodo più grande
% COSTRUZIONE DELLA MATRICE TRIDIAGONALE
h(1)=x(2)-x(1);
A(1,1)=h(1)/3;
A(1,2)=h(1)/6;
A(2,1)=h(1);
b(1)=(y(2)-y(1))/h(1)-k0;
for
i=1:n-2
h(i+1)=x(i+2)-x(i+1);
A(i+1,i+1)=2*(h(i)+h(i+1));
A(i+1,i+2)=h(i+1);
if i<n-2
A(i+2,i+1)=A(i+1,i+2);
end
b(i+1)=6*((y(i+2)-y(i+1))/h(i+1)+(y(i)-y(i+1))/h(i));
end
A(n,n-1)=h(n-1)/6;
A(n,n)=h(n-1)/3;
b(n)=k1-(y(n)-y(n-1))/h(n-1);
[L,U,risp]=lutrid(A,n);
% Fattorizzazione
LU per matrici tridiagonali
if risp==0
% RISOLUZIONE DEL SISTEMA L * U * m = b
z(1)=b(1);
% RISOLUZIONE DEL SISTEMA L * z = b
for i=2:n
z(i)=b(i)-L(i,i-1)*z(i-1);
end
m(n)=z(n)/U(n,n);
% RISOLUZIONE DEL SISTEMA U * m = z
for
i=n-1:-1:1
m(i)=(z(i)-U(i,i+1)*m(i+1))/U(i,i);
end
for
k=1:m1
if xi(k)<=x(2)
& xi(k)>=x(1)
b=xi(k)-x(1);
s(k)=y(1)+k0*b+(m(1)*(b^2))/2+(((m(2)-m(1))/h(1))*(b^3))/6;
end
for i=2:n-1
if
xi(k)<=x(i+1) & xi(k)>=x(i)
c=(y(i+1)-y(i))/h(i)-m(i)*(h(i)/3)-m(i+1)*(h(i)/6);
b=xi(k)-x(i);
s(k)=y(i)+c*b+(m(i)*(b^2))/2+(((m(i+1)-m(i))/h(i))*(b^3))/6;
end
end
end
else
s=0;
end
s=s';
if risp==0
% Algoritmo seguente utilizzato
asc=x(1):.1:x(n);
%
unicamente per tracciare il grafico
asc=asc';
[m2,h1]=size(asc);
for k=1:m2
if asc(k)<=x(2)
& asc(k)>=x(1)
b=asc(k)-x(1);
ord(k)=y(1)+k0*b+(m(1)*(b^2))/2+(((m(2)-m(1))/h(1))*(b^3))/6;
end
for i=2:n-1
if
asc(k)<=x(i+1) & asc(k)>=x(i)
c=(y(i+1)-y(i))/h(i)-m(i)*(h(i)/3)-m(i+1)*(h(i)/6);
b=asc(k)-x(i);
ord(k)=y(i)+c*b+(m(i)*(b^2))/2+(((m(i+1)-m(i))/h(i))*(b^3))/6;
end
end
end
ord=ord';
plot(asc,ord,'w',xi,s,'+r',x,y,'og');
title('Spline
cubica completa');
else
disp('il
grafico non é stato possibile disegnarlo');
end
% SPLINE CUBICA NATURALE INTERPOLANTE
function
[s,risp]=splcub(x,y,xi,n,m1);
% s
vettore soluzione dei punti da interpolare
% risp
controlla se il metodo é applicabile
% x
vettore dei nodi interpolanti
% y
valori della funzione nei nodi x
% xi
vettore dei punti in cui calcolare la funzione interpolante
% n
dimensione di x , y
% m1
dimensione di xi
% COSTRUZIONE DELLA MATRICE TRIDIAGONALE
h(1)=x(2)-x(1);
h(2)=x(3)-x(2);
A(1,1)=2*(h(2)+h(1));
A(1,2)=h(2);
A(2,1)=A(1,2);
b(1)=6*((y(3)-y(2))/h(2)+(y(1)-y(2))/h(1));
for
i=2:n-3
h(i+1)=x(i+2)-x(i+1);
A(i,i)=2*(h(i)+h(i+1));
A(i,i+1)=h(i+1);
A(i+1,i)=A(i,i+1);
b(i)=6*((y(i+2)-y(i+1))/h(i+1)+(y(i)-y(i+1))/h(i));
end
h(n-1)=x(n)-x(n-1);
A(n-2,n-2)=2*(h(n-1)+h(n-2));
b(n-2)=6*((y(n)-y(n-1))/h(n-1)+(y(n-2)-y(n-1))/h(n-2));
[L,U,risp]=lutrid(A,n-2);
% Fattorizzazione LU per matrici
tridiagonali
if risp==0
% RISOLUZIONE DEL SISTEMA L * U * m = b
z(1)=b(1);
% RISOLUZIONE DEL SISTEMA L * z = b
for
i=2:n-2
z(i)=b(i)-L(i,i-1)*z(i-1);
end
m(n-2)=z(n-2)/U(n-2,n-2);
% RISOLUZIONE DEL SISTEMA U * m = z
for
i=n-3:-1:1
m(i)=(z(i)-U(i,i+1)*m(i+1))/U(i,i);
end
for
k=1:m1
if xi(k)<=x(2)
& xi(k)>=x(1)
c=(y(2)-y(1))/h(1)-(h(1)*m(1))/6;
s(k)=y(1)+c*(xi(k)-x(1))+((m(1)/h(1))*(xi(k)-x(1))^3)/6;
end
for i=2:n-2
if
xi(k)<=x(i+1) & xi(k)>=x(i)
c=(y(i+1)-y(i))/h(i)-m(i-1)*(h(i)/3)-m(i)*(h(i)/6);
b=xi(k)-x(i);
s(k)=y(i)+c*b+(m(i-1)*(b^2))/2+(((m(i)-m(i-1))/h(i))*(b^3))/6;
end
end
if xi(k)<=x(n)
& xi(k)>=x(n-1)
c=(y(n)-y(n-1))/h(n-1)-m(n-2)*(h(n-1)/3);
s(k)=y(n-1)+c*(xi(k)-x(n-1))+(m(n-2)*(xi(k)-x(n-1))^2)/2-((m(n-2)/h(n-1))*(xi(k)-x(n-1))^3)/6;
end
end
else
s=0;
end
s=s';
if risp==0
% Algoritmo seguente utilizzato
asc=x(1):.1:x(n);
%
unicamente per tracciare il grafico
asc=asc';
[m2,h1]=size(asc);
for k=1:m2
if asc(k)<=x(2)
& asc(k)>=x(1)
c=(y(2)-y(1))/h(1)-(h(1)*m(1))/6;
ord(k)=y(1)+c*(asc(k)-x(1))+((m(1)/h(1))*(asc(k)-x(1))^3)/6;
end
for i=2:n-2
if
asc(k)<=x(i+1) & asc(k)>=x(i)
c=(y(i+1)-y(i))/h(i)-m(i-1)*(h(i)/3)-m(i)*(h(i)/6);
b=asc(k)-x(i);
ord(k)=y(i)+c*b+(m(i-1)*(b^2))/2+(((m(i)-m(i-1))/h(i))*(b^3))/6;
end
end
if asc(k)<=x(n)
& asc(k)>=x(n-1)
c=(y(n)-y(n-1))/h(n-1)-m(n-2)*(h(n-1)/3);
ord(k)=y(n-1)+c*(asc(k)-x(n-1))+(m(n-2)*(asc(k)-x(n-1))^2)/2-((m(n-2)/h(n-1))*(asc(k)-x(n-1))^3)/6;
end
end
ord=ord';
plot(asc,ord,'w',xi,s,'+r',x,y,'og');
title('Spline
cubica naturale');
else
disp('il
grafico non é stato possibile disegnarlo');
end
% CALCOLO DELLA SPLINE LINEARE INTERPOLANTE
function
s=spllin(x,y,xi,n,m);
% s vettore
soluzione dei punti da interpolare
% x vettore
dei nodi interpolanti
% y valori
della funzione nei nodi x
% xi vettore
dei punti in cui calcolare la funzione interpolante
% n dimensione
di x , y
% m dimensione
di xi
for
k=1:m
s(k)=0;
if xi(k)<=x(2)
& xi(k)>=x(1)
s(k)=s(k)+((xi(k)-x(2))/(x(1)-x(2)))*y(1);
end
for i=2:n-1
if
xi(k)>=x(i-1) & xi(k)<x(i)
s(k)=s(k)+((xi(k)-x(i-1))/(x(i)-x(i-1)))*y(i);
end
if xi(k)>=x(i)
& xi(k)<=x(i+1)
s(k)=s(k)+((xi(k)-x(i+1))/(x(i)-x(i+1)))*y(i);
end
end
if xi(k)>=x(n-1)
& xi(k)<=x(n)
s(k)=s(k)+((xi(k)-x(n-1))/(x(n)-x(n-1)))*y(n);
end
end
s=s';
plot(x,y,'w',xi,s,'+r',x,y,'og');
title('Spline
lineare passante per i nodi immessi');
% REGOLA DI HORNER PER CALCOLARE IL VALORE DEL POLINOMIO IN
UN PUNTO
% per il calcolo del valore del polinomio in un punto
conoscendo i suoi coefficienti
function
p=valpol(x,a);
% p valore
del polinomio nel punto x
% a coefficienti
del polinomio
[n,m]=size(a);
if
n==1
n=m;
end
p=0;
for
i=n:-1:1
p=p*x+a(i);
end
% CALCOLO DEGLI ZERI DI UN POLINOMIO
function [a,risp,dom]=zerpol(v1,num,err);
% a
vettore degli zeri del polinomio
% risp
controlla se il metodo é applicabile
% dom
tiene conto del fatto che si é risolto una equazione di secondo grado
% v1
coefficienti del polinomio
% num
numero di zeri reali del polinomio
% err
é l'errore massimo che si deve commettere
risp=0;
dom=0;
v3=v1;
[m,n]=size(v3);
n=m;
j1=0;
while
j1<m
if v3(1)==0
v3=v3(2:n);
% ricerca di eventuali zeri nulli
a(j1+1)=0;
n=n-1;
j1=j1+1;
else
m=j1;
end
end
[n,i]=size(v3);
if n==2
a(j1+1)=-v3(1)/v3(2);
% caso di un solo zero non nullo
end
if
n==3
Delta=sqrt((v3(2))^2-4*v3(3)*v3(1));
% caso di due zeri non nulli :
a(j1+1)=(-v3(2)-Delta)/(2*v3(3));
% solo in questo caso vengono
a(j1+2)=(-v3(2)-Delta)/(2*v3(3));
% calcolati anche gli zeri
dom=1;
%
complessi
end
if
n>3
a(j1+1)=0;
% fa in modo che non si considerino
w=zeros(num,1);
%
gli zeri ai quali é stata applicata
% la regola di Ruffini
for i=2:n
% calcolo di una maggiorazione
a(j1+1)=a(j1+1)+abs(v3(i)/v3(1));
% dello zero di modulo massimo
end
a(j1+1)=max(1,a(j1+1));
for i=1:n-1
% calcolo dei coefficienti della
v2(i)=i*v3(i+1);
% derivata prima della funzione
end
c=valpol(a(j1+1),v3);
for
j=j1:num-1
if c~=0
d=valpol(a(j+1),v2);
% prima
iterata
e=0;
for i=j1+1:j
if
w(i)==0
e=e+c/(a(j+1)-a(i));
end
end
e=d-e;
if e==0
risp=1;
else
l=a(j+1);
a(j+1)=a(j+1)-c/e;
f=valpol(a(j+1),v3);
de=valpol(a(j+1),v2);
end
while (c*f>0)
& (de*d>0) & (risp==0)
l=a(j+1);
c=f;
d=de;
% metodo di Newton a doppio passo
e=0;
for
i=j1+1:j
if w(i)==0
e=e+c/(a(j+1)-a(i));
end
end
e=d-e;
if e==0
risp=1;
else
a(j+1)=a(j+1)-2*c/e;
f=valpol(a(j+1),v3);
de=valpol(a(j+1),v2);
end
end
if j~=num-1
& risp==0
a(j+2)=a(j+1);
end
g=f;
if f~=0
& risp==0
while
(((abs(f)>err) | (abs(l-a(j+1))>err))) & (risp==0)
c=f;
l=a(j+1);
d=valpol(a(j+1),v2);
% metodo di Newton -Rapson
e=0;
for i=j1+1:j
if
w(i)==0
e=e+c/(a(j+1)-a(i));
end
end
e=d-e;
if e==0
risp=1;
else
a(j+1)=a(j+1)-c/e;
f=valpol(a(j+1),v3);
end
end
end
c=g;
else
if j~=num-1
v4(n-1)=v3(n);
for
i=n-2:-1:1
v4(i)=v3(i+1)+a(j+1)*v4(i+1);
end
v3=v4'
v4=0;
% regola di Ruffini
n=n-1;
% quando si conosce
c=valpol(a(j+1),v3);
% con precisione lo
a(j+2)=a(j+1);
% zero del polinomio
w(j+1)=1;
v2=0;
for
i=1:n-1
v2(i)=i*v3(i+1);
end
end
end
end
end
a=a';
Next: Risoluzione di sistemi
lineari e fattorizzazioni di matrici (pag.1) Up:Introduzione
Previous: Function
da 21 a 30
Intro Gen: Introduzione
Generale Home: Home
page
Vito Marinelli
8-5-2000