Rieccomi,
ecco la soluzione, con un codice eseguibile:
Come scritto prima l'equazione discretizzata è, con la convenzione $u_i \approx u(x_i)$:
$ \frac{u_{i+1} - 2*u_i + u_{i-1}}{h^2} + \sin(x_i)=0 $ (***)
Come appare evidente dalla precedente scrittura, il problema si riduce alla risoluzione di un sistema lineare di dimensione $m$, dove $m$ è il numero di nodi utilizzati. Nota che è
lineare, e non serve il metodo di Newton, come sicuramente avrai visto a lezione. Quindi, se chiamiamo $D_2$ l'operatore differenziale che discretizza la derivata seconda della soluzione, si ha che il problema si riscrive come:
$F(\vec{u})=D_2 \vec{u} + \sin(x)=\vec{0}$
Nel nostro caso, $D_2$ è
rappresentato semplicemente da una matrice di dimensione $m$,
tridiagonale. Pertanto, si tratta di trovare la soluzione del sistema lineare $D_2 \vec{u} = -\sin(x)$
Nodo $x=1$Si deve imporre che $u_1=0$. E' sufficiente cambiare la prima riga, ricordando che per il primo nodo ($x_1=0$), si ha automaticamente $\sin(x_1)=\sin(0)=0$.
Nodo $x=2 \pi$ Questo è nel messaggio precedente. In particolare, si ricava $u_{m+1}=u_{m-1}+2h$.
Sostituendo questa espressione in (***) si ha:
$\frac{2}{h^2} u_{m-1} - \frac{2}{h^2}u_m + \frac{2}{h}=0$
e questa sarà l'ulima riga del tuo sistema da risolvere.
A questo punto devi risolvere il ssitema lineare $D_2 \vec{u} = - \sin(x)$, dove $x$ ovviamente è un vettore, e dunque, $\sin(x)$ pure è un vettore di $m$ componenti. Ricorda che sia $D_2$ che $b$ sono stati modificati nelle righe corrispondenti alle condizioni ai bordi.
ImplementazioneQuesto è il codice. Ho mostrato un grafico di convergenza, dove si vede che l'ordine è effettivamente due. Inoltre, nota come abbia confrontato ogni soluzione ottenuta con $m$ passi con la soluzione analitica, che è ovviamente $u(x)=\sin(x)$.
- Codice:
clear all
close all
counter=0;
mrange=2.^(4:10)+1;
for m=mrange
counter=counter+1;
h=2*pi/(m-1);
x=linspace(0,2*pi,m)';
D2 = toeplitz(sparse([1,2],[1,1],[-2,1]/(h^2),m,1));
b=sin(x);
%condizioni ai bordi
%nodo x=1
D2(1,1:2)=[1,0];
%nodo x=2pi
D2(m,m-1:m)=[2/h^2,-2/h^2]; %termine da cambiare che coinvolge l'incognita
b(m)=2/h; %cambio ultima riga del termine noto
u=D2\(-b); %risoluzione sistema lineare
err(counter)=norm(u-sin(x),Inf); %calcolo l'errore misurato rispetto alla soluzione analitica
end
loglog(mrange,err,'*',mrange,err(end)*(mrange/mrange(end)).^(-2),'g')
legend('Errore','Ordine 2')
title('Error plot')
Il grafico prodotto è il seguente
Spero ti sia stato utile, ciao!