Ordine di convergenza della spline cubica, esercizio

Angus1956
Consideriamo il seguente problema:




il mio codice è questo:

Script:

m=10;
a=-pi/2; b=pi/2;
g=@(t)(t.*cos(t).*sin(t));
f = @(x)( (x(:,end)+pi/2)/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,10000);
for tt=t
    k=k+1; 
    x(k,:)=linspace(-pi/2,tt,2*m+1);
end
figure(1)
plot(t,f(x),'r')
hold on
pause
t=linspace(a,b,10000)';
for n=[4 6 8]
     x1=linspace(a,b,n+1);
     xcap=cos( (2*(0:n)+1)./(2*n+2)*pi)';
     x2=(a+b)/2+(b-a)/2*xcap;
     for k=1:n+1
        xf1=linspace(-pi/2,x1(k),2*m+1);
        xf2=linspace(-pi/2,x2(k),2*m+1); 
        y1(k)=f(xf1);
        y2(k)=f(xf2);
     end
    s=spline(x1(:),[0;y1(:);0],t);
    c = get_polyn(x2,y2');
    yp= polyval(c,t);
    figure(1)
    plot(t,s,'b')
    plot(t,yp,'g');
    title('Grafico errore')
    legend('f','spline cubica','polinomio interpolante nei nodi di Chebyshev')
    xlabel('asse x')
    ylabel('asse y')
    pause
    figure(2)
    plot(t,f(x)-yp);
    title('Funzione Errore di interpolazione')
    xlabel('valutazioni in t')
    ylabel('Errore con nodi di Chebyshev')
    pause
end


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;


Volevo sapere come fare a trovare in questo caso l'ordine di convergenza $p$ (che di solito è $p=(1/log(1/2))*log(e_k)/log(e_(k-1))$, nel caso in cui si abbia che $h=(b-a)/2^k$):Poi volevo chiedere se l'errore del polinomio interpolante nei nodi di Chebyshev andasse bene, mi esce questo grafico (in blu $n=4$, in arancio $n=6$ e in giallo $n=8$):


Risposte
Angus1956
Io pensavo di modificare così (per trovare l'ordine di convergenza $p$):
for n=[4 6 8]
x1=(a:(b-a)/2^log2(n):b)';
...
 Er=max(abs(f(t)-s));
 if n~=4
 p=log2(Er/Er_old);
 disp([n,p]);
 end
Er_old=Er;
end

Martino
[xdom="Martino"]A giudicare dal testo, questo esercizio sarà valutato e riceverà un voto. Per cui chiudo.[/xdom]

Rispondi
Per rispondere a questa discussione devi prima effettuare il login.