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...$)

[feddy]

Dal tuo codice noto che hai scritto $x^{\star}=0.5$, mentre $x^{\star}$ va determinato (...in modo tale che $F_{x^*}(x)$ etc)

Quello che intendo dire è che mi sembra tu abbia interpolato direttamente in $0.5$ invece che in $x^{\star}$

[Angus1956]

"feddy":

Quello che intendo dire è che mi sembra tu abbia interpolato direttamente in $0.5$ invece che in $x^{\star}$

No, non ho interpolato in $0.5$, infatti non ho scritto $y(9)=f(0.5)$ ma ho scritto $y(9)=f(0.5)pmsqrt(10^-3)$ ovvero imposto che in $0.5$ il polinomio deve rispettare la condizione imposta da $F$, cioè mi son trovato il valore di $p$ in $0.5$ con la condizione che $F$ vale $0.001$ in $0.5$ (non so se hai capito), altrimenti non saprei bene come fare.

[Angus1956]

C'è so solo che $F$ si annulla dove $p$ interpola $f$ e so quanto vale in $0.5$, per trovare $x^**$ un idea sarebbe quella di vedere gli zeri di $F$ ma come faccio a beccare proprio $x^**$?

[feddy]

La traccia vuole farti usare la stima dell'errore di interpolazione

[Angus1956]

"feddy":La traccia vuole farti usare la stima dell'errore di interpolazione

$||f^((n+1))(x)||_{infty}/((n+1)!)||omega(x)||_infty $ dove $omega(x)=(x-x_0)*...*(x-x_7)*(x-x^**)$? Scusa a cosa mi servirebbe?

[feddy]

La funzione $F$ è proprio l'errore (o meglio, il suo quadrato) per cui vale quella stima che hai scritto. Usando quella stima puoi trovare un'equazione con incognita $x^{\star}$

[Angus1956]

scusa non è tropo grossolana come approssimazione ( e quindi $x^**$ non è precisa)? Comunque quindi mi basta calcolare l'errore al quadrato in $0.5$ e trovare $x^**$, ma la derivata decima la devo calcolare a mano (sarebbe tipo $-512 (-5 cos(2 x) + x sin(2 x))$ se non sbaglio)?

[feddy]

Grossolana o no va deciso in base all'accuratezza che hai raggiunto con i punti precedenti. Per quanto riguarda la derivata decima, in un modo o nell'altro la devi calcolare, o stimare! Infatti, la tua formula dell'errore la derivata è sbagliata perchè la derivata $n+1$-esima non è valutata in $x$. Controlla la teoria.

[Angus1956]

"feddy":Grossolana o no va deciso in base all'accuratezza che hai raggiunto con i punti precedenti. Per quanto riguarda la derivata decima, in un modo o nell'altro la devi calcolare, o stimare! Infatti, la tua formula dell'errore la derivata è sbagliata perchè la derivata $n+1$-esima non è valutata in $x$. Controlla la teoria.

Si scusa intendevo in un certo $\xi \in[a,b]$. Comunque ho visto un commento (riguardo l'esercizio) della professoressa e c'è scritto questo:
"
Possibile handle per richiamare il metodo di bisezione, per es.

FF=@(t)( (f(0.5) - polyval(polyfit([x,t],f([x,t]),8),0.5) )^2-0.001);

(qui viene proposto polyfit, ma si può fare con get_polin)"

Quindi si trova $x^**$ come radice di $F$ (probabilmente vedendo dal grafico di $F$ dove si trova $x^**$ e usa la bisezione su un intervallo che contiene solo $x^**$ come radice di $F$).

[Angus1956]

Però non ho capito perchè usa come polinomio

polyfit([x,t],f([x,t]),8)

in teoria questo è il polinomio interpolante $f$ nei nodi $x_0,...,x_7$ mentre nella scrittura di $F$ dovrebbe essere quello che interpola anche $x^**$? (inoltre se non ricordo male la sintassi di polyfit è del tipo $polyfit(x,y,n)$ dove $x$ è il vettore dei nodi in cui interpolare $f$, poi $y$ sono i valori di $f$ nei nodi componenti di $x$ e infine $n$ è il grado del polinomio interpolante, ma allora perchè come nodi ha usato $[x,t]$, forse per la dipendenza di $F$ da $t$?Oltre il fatto che $[x,t]$ mi da un errore di concatenazione dato che $x$ è un vettore e $t$ uno scalare...)

[feddy]

Sì, vuole trovare lo zero di quella funzione, e per farlo devi usare un metodo opportuno. Non conosco bene la sintassi MatLab, dovrei provare sinceramente. Quel termine è l'interpolante nei nodi $x_0$ fino a $x^{\star}$, che nella function è scritto come $t$ in quanto incognita del problema.

[Angus1956]

"feddy": Quel termine è l'interpolante nei nodi $x_0$ fino a $x^{\star}$, che nella function è scritto come $t$ in quanto incognita del problema.

a ok capito, ma allora si dovrebbe scrivere come $[x]$

[feddy]

Il problema dovrebbe essere che x è un vettore colonna, al quale vuoi aggiungere t, e questo non è possibile. Potresti aggiungere un apice ad x all'interno della funzione.

[Angus1956]

"feddy":Il problema dovrebbe essere che x è un vettore colonna, al quale vuoi aggiungere t, e questo non è possibile. Potresti aggiungere un apice ad x all'interno della funzione.

In teoria scusa ma se io ho già un vettore e poi uno scalare se voglio ottenere il vettore dato dal vettore di partenza e l'ultima componente lo scalare non basta fare $[x]$?

[feddy]

Quando hai questo tipo di dubbi ti conviene provare direttamente tu stesso sulla command window:

>> x = linspace(0,1,5)

x =

         0    0.2500    0.5000    0.7500    1.0000

>> t = 44;
>> [x,t]

ans =

         0    0.2500    0.5000    0.7500    1.0000   44.0000

Come vedi il valore di $t$ è stato appeso alla fine di $x$, mentre se fai

>> [x]

ans =

  Columns 1 through 11

     0     1     2     3     4     5     6     7     8     9    10

  Columns 12 through 22

    11    12    13    14    15    16    17    18    19    20    21

  Columns 23 through 33

    22    23    24    25    26    27    28    29    30    31    32

  Columns 34 through 44

    33    34    35    36    37    38    39    40    41    42    43

  Column 45

    44

hai scritto un range equispaziato (come quando scrivi un for-loop) da 0 a 44.

[feddy]

Ho cambiato un po' il tuo codice, usando il suggerimento della prof.

FF

è la funzione di cui vuoi trovare lo zero $x^{\star}$. Ho prima plottato un po' di valori per scegliere un guess iniziale, e ho usato la funzione di MatLab

fzero

per calcolarne lo zero. Dopodichè, ho usato quanto avevi già fatto tu.

clear all
close all

f=@(x)(x.*sin(x).*cos(x));
alpha=0;
beta=2*pi;
figure(1)
fplot(f,[alpha,beta],'r--')
hold on

x=linspace(alpha,beta,8)';
plot(x, f(x),'*');
hold on
y = f(x);
a = get_polyn(x,y);
t = linspace(alpha,beta,100);
yp=polyval(a,t);
plot(t,yp,'o');

%Punto 2
FF=@(xstar)  (f(0.5) - polyval(polyfit([x',xstar],f([x',xstar]),8),0.5) ).^2 - 1e-3;

xstarrange=linspace(alpha,beta,50);
ystarrange= zeros(1,numel(xstarrange));
for i=1:numel(xstarrange)
    ystarrange(1,i)=FF(xstarrange(i));
end
% plot(xstarrange, ystarrange,'o');
xstar=fzero(FF,0.1); %trova lo zero, guess iniziale x0 = 0.1

x(9)=xstar(end);
y(9)=f(x(9));
a = get_polyn(x,y);
t = linspace(alpha,beta,100);
yp=polyval(a,t);
hold on
plot(t,yp,'k--');
title('Grafico e polinomio interpolante')
legend('f','nodi','polinomio interpolante','polinomio interpolante con x^*')
xlabel('asse x')
ylabel('asse y')


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;
end

[Angus1956]

invece di plottare i valori e usare fzero, non si può usare tipo il metodo di bisezione e ti fai i grafico di $F$ per capire dove sta $x^**$ (vai a esclusione con gli altri nodi) e quindi poi ti trovi l'intervallo dove applicare la bisezione?

[feddy]

La ricerca dello zero (con bisezione ad esempio) va applicata alla funzione FF, non F come hai scritto nel messaggio qui sopra. Ho plottato i valori di FF per comodità sinceramente, è solamente un'opzione per trovare un intervallo opportuno.

Scusa ma non capisco cosa intendi sinceramente dicendo:

"andreadel1988":non puoi usare tipo il metodo di bisezione e ti fai i grafico di F per capire dove sta x∗ (vai a esclusione con gli altri nodi) e quindi poi ti trovi l'intervallo dove applicare la bisezione?

è come se avessi scritto "non puoi usare il metodo di bisezione [...] per trovare l'intervallo dove applicare la bisezione?". Intendi dire di vedere per quali coppie di estremi il prodotto è negativo e dunque usare bisezione? Un metodo vale l'altro, una volta che hai trovato un intervallo dove cambia il segno va bene.

[Angus1956]

"feddy":La ricerca dello zero (con bisezione ad esempio) va applicata alla funzione FF, non F come hai scritto nel messaggio qui sopra. Ho plottato i valori di FF per comodità sinceramente, è solamente un'opzione per trovare un intervallo opportuno.

Scusa ma non capisco cosa intendi sinceramente dicendo:
[quote="andreadel1988"]non puoi usare tipo il metodo di bisezione e ti fai i grafico di F per capire dove sta x∗ (vai a esclusione con gli altri nodi) e quindi poi ti trovi l'intervallo dove applicare la bisezione?

è come se avessi scritto "non puoi usare il metodo di bisezione [...] per trovare l'intervallo dove applicare la bisezione?". Intendi dire di vedere per quali coppie di estremi il prodotto è negativo e dunque usare bisezione? Un metodo vale l'altro, una volta che hai trovato un intervallo dove cambia il segno va bene.[/quote]
Scusa mi sono espresso male, intendevo dire che prima facevo il grafico di $FF$ per trovare un intervallo opportuno dove applicare il metodo di bisezione (ovvero il prodotto delle due ascisse calcolate in $f$ fa $0$ e che non sia un nodo tra $x_0,...x_7$ che sono dati) e poi applica il metodo di bisezione in questo intervallo. Una domanda, quanto dovrebbe venire $x^**$?

[Angus1956]

Comunque a me il grafico di $FF$ viene questo:



,ho usato questo codice:

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;
figure(2)
fplot(FF,[alpha,beta],'b')