Polinomio composito e area di una funzione integrale

[Angus1956]

Consideriamo il seguente esercizio:



Ho io scritto questo codice

Script:

g=@(t)(t.*(sin(t)).^2.*exp(-t));
m=10;
a=-2;
b=1;
F=@(x)(x-(x(:,end)-0)/6/m.*( g(x(:,1))+2*sum(g(x(:,3:2:2*m)),2)+4*sum(g(x(:,2:2:2*m)),2)+g(x(:,2*m+1))) );
k=0;
t=linspace(a,b,100);
for tt=t
    k=k+1; 
    x(k,:)=linspace(0,tt,2*m+1);
end
figure(1)
hold on
plot(t,F(x),'k')
pause
n=3;
A=0;
x = linspace(a,b,m+1);
for j=1:m
    x_in=linspace(x(j),x(j+1),n+1)'; 
    for i=1:length(x_in)
        y(i) = x_in(i)-Cavalieri_Simpson(-1,x_in(i),m,g);
    end
    a = get_polyn(x_in,y');
    t = linspace(x(j),x(j+1),100);
    yp= polyval(a,t);
    plot(t,yp,'r');
    A=A+Cavalieri_Simpson_polinomi(x(j),x(j+1),m,a);
    pause
end
title('Grafico')
legend('F','polinomio composito')
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 Cavalieri_Simpson_polinomi:

function I2m=Cavalieri_Simpson_polinomi(a,b,m,f)
H=(b-a)/m;

x=linspace(a,b,2*m+1);

I2m=H/6*( polyval(f,x(1))+ 2*sum(polyval(f,x(3:2:2*m))) + 4*sum(polyval(f,x(2:2:2*m))) + polyval(f,x(2*m+1)));

end

Il tutto funziona, però quando va a fare il grafico della funzione $F$ me lo da sbagliato, come posso risolvere senza modificare la base della formula di quadratura?
.

[feddy]

Cosa intendi con sbagliato?

[Angus1956]

Avevo pensato di riscrivermi $x=\int_{-1}^{x}f(t)dt$ in modo tale che $F(x)=\int_{-1}^{x}f(t)-t(sint)^2e^(-t)dt$, però non credo che l'equazione integrale $x=\int_{-1}^{x}f(t)dt$ in $f(t)$ sia facilmente risolvibile (non ne ho mai fatta una ahhaha)

[Angus1956]

"feddy":Cosa intendi con sbagliato?

Questo è il grafico di $F$ che mi esce dallo handle:



come vedi mi viene tipo un fascio di funzioni e pure sbagliato poichè in $0$ non deve valere $0$ mentre il polinomio interpolante è giusto

[feddy]

Se dopo ho un po' di tempo do un'occhiata al codice. A livello concettuale è lo stesso problema che avevi nel post precedente, ossia una funzione che all'interno contiene un integrale, ma ammetto di non aver letto il tuo codice.

[Angus1956]

"feddy":Se dopo ho un po' di tempo do un'occhiata al codice. A livello concettuale è lo stesso problema che avevi nel post precedente, ossia una funzione che all'interno contiene un integrale, ma ammetto di non aver letto il tuo codice.

Eh si, però devo per forza usare questa forma..., grazie.

[feddy]

Capisco, scusami ma non lo ricordavo e nella traccia non c'è scritto... Il problema è che la tua $F$ è definita in modo tale da prendere un linspace (per utilizzare la formula di quadratura), ma il primo addendo della funzione non deve essere un linspace, ma solamente uno scalare ($x$).

[Angus1956]

"feddy":Capisco, scusami ma non lo ricordavo e nella traccia non c'è scritto... Il problema è che la tua $F$ è definita in modo tale da prendere un linspace (per utilizzare la formula di quadratura), ma il primo addendo della funzione non deve essere un linspace, ma solamente uno scalare ($x$).

Si esatto, comunque non devi scusarti di niente, non ho capito perchè ti dovresti scusare, se parli della forma ahimè come ti ho detto non mi è concesso usare codici troppo fuori da quelli fatti nel corso (non posso definire funzioni con altre funzioni anche da me create ma solo con gli handle o valutandole nei nodi)... però l'idea di usare la roba composita è direi diretta, poi bisogna vedere come usarla bene.

[feddy]

Se vuoi usare questa forma ti è sufficiente aggiustare la $F$ in questo modo:

F=@(x)(x(end)-(x(:,end)-0)/6/m.*( g(x(:,1))+2*sum(g(x(:,3:2:2*m)),2)+4*sum(g(x(:,2:2:2*m)),2)+g(x(:,2*m+1))) );

Nota l'x(end) all'inizio. Il resto è uguale

[Angus1956]

"feddy":Se vuoi usare questa forma ti è sufficiente aggiustare la $F$ in questo modo:

F=@(x)(x(end)-(x(:,end)-0)/6/m.*( g(x(:,1))+2*sum(g(x(:,3:2:2*m)),2)+4*sum(g(x(:,2:2:2*m)),2)+g(x(:,2*m+1))) );

Nota l'x(end) all'inizio. Il resto è uguale

Se posso chiedere vai prendere tutto il vettore $x$ per omogeneizzare il fatto che sia una funzione definita vettorialmente? Devo essere sincero questa formula non l ho mai capita appieno proprio perchè utilizza vettori e tutto su una funzione e nel corso non è stato spiegato minimamento come funziona una funzione su matlab a piu variabili (ovviamente so come funzionano le funzioni a piu variabili, intendo direi cosa fa matlab in questo caso per definire la funzione ) però vorrei capire bene su matlab cosa fa, passo dopo passo.

[Angus1956]

Comuqnue ho trovato un altro errore, ovvero al posto di $0$ andavo sostituito $-1$ poichè l'integrale va da $-1$ a $x$ (mentre nel codice da dove l'avevo copiato andava da $0$ a $x$), perciò l'ho riscritto così:
Script:

g=@(t)(t.*(sin(t)).^2.*exp(-t));
m=10;
a=-2;
b=1;
F=@(x)(x(end)-(x(:,end)+1)/6/m.*( g(x(:,1))+2*sum(g(x(:,3:2:2*m)),2)+4*sum(g(x(:,2:2:2*m)),2)+g(x(:,2*m+1))) );
k=0;
t=linspace(a,b,100);
for tt=t
    k=k+1; 
    x(k,:)=linspace(-1,tt,2*m+1);
end
figure(1)
hold on
plot(t,F(x),'k')
pause
n=3;
A=0;
x = linspace(a,b,m+1);
for j=1:m
    x_in=linspace(x(j),x(j+1),n+1)'; 
    for k=1:n+1
        xF=linspace(-1,x_in(k),2*m+1); 
        y(k)=F(xF);
    end
    a = get_polyn(x_in,y');
    t = linspace(x(j),x(j+1),100);
    yp= polyval(a,t);
    plot(t,yp,'r');
    A=A+Cavalieri_Simpson_polinomi(x(j),x(j+1),m,a);
    pause
end
title('Grafico')
legend('F','polinomio composito')
xlabel('asse x')
ylabel('asse y')

però il polinomio interpolante non interpola per niente $F$ nei nodi, anche se è molto simile a $F$, questi sono i due grafici:



non capisco cosa sia sbagliato, sembra tipo traslato il polinomio interpolante rispetto a $F$...

[feddy]

Esatto, ho voluto evitare il fatto che il primo addendo risultasse una matrice. x è una matrice di taglia 100 x 21. Ad ogni riga di x fissata corrisponde un vettore coi punti necessari per applicare la formula di quadratura da $0$ ad un certo valore compreso tra $-2$ e $1$.

EDIT: Si, -1 è corretto, non avevo guardato

[Angus1956]

Ho provato a vedere il valore in $0$ di $F$ ma non rida:



invece dovrebbe ridare quaclosa del tipo $0,44439$, è un errore di approssimazione di cavalieri-simpson o è proprio sbagliato? (anche se come errore di approssimazione mi sembra molto grossolano, inoltre mi viene il dubbio che sia questo il problema per cui il polinomio non interpola $F$, quindi è il polinomio a essere esatto non $F$ aahahahhaha)

[Angus1956]

In effetti il polinomio è corretto:



è $F$ che è sbagliata

[feddy]

Usando la tua $F$ e interpolando con polyfit + polyval (quindi senza il tuo codice) il tutto funziona.

[Angus1956]

"feddy":Usando la tua $F$ e interpolando con polyfit + polyval (quindi senza il tuo codice) il tutto funziona.

A me sembra tipo che la mia $F$ sia tralsata di $1$ sopra l asse delle y ma non capisco perchè (infatti se tipo in $0$ ci tolgo 1 mi viene $0,44439$)

[Angus1956]

"feddy":Usando la tua $F$ e interpolando con polyfit + polyval (quindi senza il tuo codice) il tutto funziona.

Comunque polyval posso usarlo per valutare i polinomi, polifit invece non lo posso usare anche se so cos'è ma non ci è consentito, comunque il problema non è nel polinomio interpolante ma in $F$ per quello che compare dai grafici

[feddy]

Perdonami, la mia $F$ era sbagliata, l'avevo data per buona . L'errore è dunque nella funzione in cui interpoli. Anche nel caso di una $F$ sbagliata, avresti dovuto trovare l'interpolante. Prova a testare la tua funzione interpolante in casi più semplici per vedere dove è il problema.

[Angus1956]

"feddy":Perdonami, la mia $F$ è sbagliata e la tua è giusta. Infatti in -2 la mia $F$ ha un valore errato come puoi notare. L'errore è dunque nel modo in cui interpoli.

Scusami come fa a essere sbagliato il polinomio interpolante se esso ridà in $0$ mentre $F$ non ridà in $0$? Se provi tipo a usare un calcolatore per vedere quanto dovrebbe fare $F$ in $0$ ti viene tipo $0.44439$ che è giusto per il polinomio interpolante e sbagliato per $F$ che è tipo traslata $1$ di sopra (parlo dei miei grafici che ho mandato), ad esmepio io ho usato wolphrama per calcolare $F$ in $0$ e viene:

[feddy]

Sostituisci la tua $F$ con questa:

F=@(x)(x(:,end)-(x(:,end)+1)/6/m.*( g(x(:,1))+2*sum(g(x(:,3:2:2*m)),2)+4*sum(g(x(:,2:2:2*m)),2)+g(x(:,2*m+1))) );

Ora dovrebbe tornare il tutto. Confermi?

[Angus1956]

"feddy":Sostituisci la tua $F$ con questa:

F=@(x)(x(:,end)-(x(:,end)+1)/6/m.*( g(x(:,1))+2*sum(g(x(:,3:2:2*m)),2)+4*sum(g(x(:,2:2:2*m)),2)+g(x(:,2*m+1))) );

Ora dovrebbe tornare il tutto. Confermi?

Si è vero!!! hai messo $x( :,end)$ in pratica, quindi una matrice, mamma mia non ci avrei pensato grazie mille