Codice matLab polinomi interpolanti

[Angus1956]

Considerato questo problema sui polinomi interpolanti:





Volevo sapere se il mio codice potesse andare bene (soprattutto il punto 2):

Script

f=@(x)(x.*sin(x).*cos(x));
alpha=0;
beta=2*pi;
figure(1)
fplot(f,[alpha,beta],'r--')
hold on
pause
x=linspace(alpha,beta,8)';
plot(x, f(x),'*');
y = f(x);
a = get_polyn(x,y);
t = linspace(alpha,beta,100);
yp=polyval(a,t);
plot(t,yp,'k');
pause
x(9)=0.5;
y(9)=f(0.5)+sqrt(1e-3);
a = get_polyn(x,y);
t = linspace(alpha,beta,100);
yp=polyval(a,t);
plot(t,yp,'b');
title('Grafico e polinomio interpolante')
legend('f','nodi','polinomio interpolante','polinomio interpolante con x^*')
xlabel('asse x')
ylabel('asse y')

Funzione:

function a = get_polyn(x,y)
n=length(x)-1;
V(1:n+1,n+1)=ones(n+1,1);
for j=n:-1:1
 V(1:n+1,j) = x.*V(1:n+1,j+1);
end
a = V\y;

Grazie (ho visto anche graficamente e me lo interpola in due nodi in più ovvero $0,7483...$ e $2,5908....$ non so se sia corretto in quanto ne chiedeva uno però mi sembra tutto corretto..., ho provato anche a mettere

y(9)=f(0.5)-sqrt(1e-3);

e anche qui mi veniva due nodi in più ovvero $0,2780...$ e $2.9487...$)

[Angus1956]

Che poi se rifacessi il codice iniziale della discussione mettendo

y(9)=f(0.5)-sqrt(1e-3);

, il polinomio interpolerebbe $f$ in $0.278$ che sarebbe $x^**$, quindi poteva andare bene...

[Angus1956]

Comunque il mio codice finale è questo:

Script:

f=@(x)(x.*sin(x).*cos(x));
alpha=0;
beta=2*pi;
figure(1)
fplot(f,[alpha,beta],'r--')
hold on
pause
x=linspace(alpha,beta,8)';
plot(x, f(x),'*');
y = f(x);
a = get_polyn(x,y);
r = linspace(alpha,beta,100);
yp=polyval(a,r);
plot(r,yp,'k');
pause
FF=@(t)  (f(0.5) - polyval(polyfit([x;t],f([x;t]),8),0.5) ).^2 - 1e-3;
kmax=500;
tolla=1e-8;
tollr=0;
tollf=1e-8;
b=0.1;
c=0.3;
[ind,x(9)]=bisezione(FF,b,c,tolla,tollr,tollf,kmax);
y(9)=f(x(9));
a = get_polyn(x,y);
yp=polyval(a,r);
plot(r,yp,'b');
title('Grafico e polinomio interpolante')
legend('f','nodi','polinomio interpolante','polinomio interpolante con x^*')
xlabel('asse x')
ylabel('asse y')

Funzione get_polyn:

function a = get_polyn(x,y)
n=length(x)-1;
V(1:n+1,n+1)=ones(n+1,1);
for j=n:-1:1
 V(1:n+1,j) = x.*V(1:n+1,j+1);
end
a = V\y;

Funzione bisezione:

function [ind,x,y,k,afinal,bfinal,k_1,res]=bisezione(f,a,b,tolla,tollr,tollf,kmax)
f_a=f(a);
ind=-1;
for k=1:kmax
    if k>1
        x_0=x;
    end
    x=(a+b)/2;
    if x>=1e-9
        k_1=k;
    end
    y=f(x);
    res(k)=abs(y);
    % criterio di arresto su f
    if abs(y)<=tollf, ind=2; return, end
    % criterio di arresto su [a,b]
    if b-a<= tolla+tollr*abs(a), ind=1; return, end
    if f_a*y<0
        b=x;
    else   
        a=x;
        f_a=y;
    end
    afinal=a;
    bfinal=b;
end

l'unica cosa che non sono riuscito a fare per bene e a fare il grafico di $FF$ per trovare l'intervallo giusto della bisezione che ho fatto a occhio vedendo il guest che hai dato tu, però boh non mi esce bene il grafico di $FF$...

[feddy]

Puoi tirare una linea orizzontale sul plot per $y=0$ e zoommare per vedere dove cambia segno. Io prima, una volta ottenute le y corrispondenti a FF avevo guardato min(y) e max(y). Siccome uno dei due era negativo, ho guardato nel plot dove erano circa quei valori e ho scelto il guess. Quando hai numeri molto piccoli bisognerebbe utilizzare un plot logaritmico. Il problema è che per valori negativi non è definito, dunque non riesci davvero a visualizzare immediatamente, ma ci sono altri trucchetti in questo caso.

[Angus1956]

"feddy":Puoi tirare una linea orizzontale sul plot per $y=0$ e zoommare per vedere dove cambia segno. Io prima, una volta ottenute le y corrispondenti a FF avevo guardato min(y) e max(y). Siccome uno dei due era negativo, ho guardato nel plot dove erano circa quei valori e ho scelto il guess. Quando hai numeri molto piccoli bisognerebbe utilizzare un plot logaritmico. Il problema è che per valori negativi non è definito, dunque non riesci davvero a visualizzare immediatamente, ma ci sono altri trucchetti in questo caso.

Più che altro non capisco se la mia $FF$ sia giusta graficamente, poi per la storia di come trovare dove fa $0$ di solito vedo i punti dove fa negativo e poi positivo.

[feddy]

La FF è molto simile alla funzione errore, o meglio: è la funzione errore traslata di 1e-3. Per vedere che sia giusta, o perlomeno ragionevole, quello che deve succedere è che in corrispondenza dei nodi di interpolazione sia molto piccola (traslata di 1e-3), ed è quello che succede.

[Angus1956]

"feddy":La FF è molto simile alla funzione errore, o meglio: è la funzione errore traslata di 1e-3. Per vedere che sia giusta, o perlomeno ragionevole, quello che deve succedere è che in corrispondenza dei nodi di interpolazione sia molto piccola (traslata di 1e-3), ed è quello che succede.

Ok, grazie, io però ancora precisamente capito perchè ha definito questa funzione $FF$ che non corrisponde con $F_{x^**}$ giusto? C'e io avevo pensato che si trovasse lo zero di $F_{x^**}$ però la funzione $FF$ che ha definito non mi sembra sia la $F_{x^**}$ definita nel problema o forse non la vedo io che si scrive così...

[feddy]

Se vuoi risolvere $f(x)=a$, trasformi il problema nella ricerca dello zero una nuova funzione $ff(x)$ definita come $$ff(x):=f(x)-a$$

La logica è la stessa, solo che qui l'incognita è più complessa ed è racchiusa nella valutazione dell'interpolante.

[Angus1956]

"feddy":Se vuoi risolvere $f(x)=a$, trasformi il problema nella ricerca dello zero una nuova funzione $ff(x)$ definita come $$ff(x):=f(x)-a$$

La logica è la stessa, solo che qui l'incognita è più complessa ed è racchiusa nella valutazione dell'interpolante.

A ok quindi lei vuole risolvere $F_{x^**}(0.5)=0.001$, io pensavo volesse risolvere $F_{x^**}(x^**)=0$ c'è trovarsi da subito $x^**$ come zero di $F_{x^**}$

[Angus1956]

"andreadel1988":
A ok quindi lei vuole risolvere $F_{x^**}(0.5)=0.001$, io pensavo volesse risolvere $F_{x^**}(x^**)=0$ c'è trovarsi da subito $x^**$ come zero di $F_{x^**}$

Quindi allora si è scritto $p(0.5)$ come $polyval(polyfit([x;t],f([x;t]),8),0.5) $ dove appunto ha messo la dipendenza dall'ultimo nodo, chiamandolo $t$ poichè non lo conosceva e quest'ultimo se lo è trovato come zero della funzione $0=F_{x^**}(0.5)-0.001=(f(0.5)-p(0.5))-0.001$. Un piccolo chiarimento: ha scritto $polyfit([x;t],f([x;t]),8)$ invece che $polyfit([x;t],f([x;t]),9)$ poichè in questo caso il nono nodo non lo abbiamo, ovvero polyfit invece di ridare i coefficenti del polinomio in questo caso rida una funzione in $t$ (?

[feddy]

Ha messo $8$ e non $9$ perchè vuole valutare proprio il polinomio determinato prima, che usa appunto $8$. Come hai scritto, il nono nodo non lo abbiamo

[Angus1956]

"feddy":Ha messo $8$ e non $9$ perchè vuole valutare proprio il polinomio determinato prima, che usa appunto $8$. Come hai scritto, il nono nodo non lo abbiamo

Si si è che mi stavo rifacendo alla sintassi di polyfit, che vuole in entrata un certo numero di nodi dati dalle componenti di un vettore $x$ e come $n$ vuole lo stesso numero di componenti del vettore $x$, solo che in questo caso avevamo $[x;t]$ un vettore di $9$ componenti e invece $n=8$, ma in questo caso perchè definiamo una funzione in funzione del nono nodo (che dobbiamo determinare) e quindi calcola il polinomio nei restanti $8$ nodi e $t$ lo usa come variabile (credo).