Ciao, benvenuto/a sul forum. Consiglio, come da
Regolamento, di scrivere in LaTeX. Vedi come si scrivono le formule:
formuleil problema è lineare, quindi si potrebbe anche usare la soluzione analitica per confrontare l'errore, ma veniamo al dunque:
con $u_i$ si intende $u (x_i)$, dove $x_i$ è il nodo i-esimo. Mentre con $h=\frac{b-a}{m}$ si intende il passo di discretizzazione (costante in questo caso), mentre $m$ è il numero dei nodi.
Condizione di NeumannPer le condizioni al bordo in questo caso la scelta è usare un
nodo fantasma:
click e
clack (il primo è in inglese, all'inizio può sembrare un po' complicato, il secondo è in italiano [vai a pag.6]).
In poche parole, questa strategia consiste nell'usare lo schema alle differenze finite per approssimare la derivata prima
centrandole nel bordo $x_1=0$. Ossia:
$\frac{u(x_2) - u(x_0)}{2h}=u'(x_1)$
Questa approssimazione è del
secondo ordine.
In spoiler un'altra strategia
Testo nascosto, fai click qui per vederlo
Se avessi voluto, avresti potuto utilizzare differenze finite
non centrate. L'idea è semplice: non si introduce alcun nodo fantasma, ma si considerano i nodi a destra o a sinistra, a seconda del caso.
In questo problema, al nodo $x_1=0$, avremo
$\frac{u(x_2) - u(x_1)}{h}=u'(x_1)$
Purtroppo è un'approssimazione del primo ordine.
Se avessimo voluto quella del secondo ordine, si ha
$u'(x_1)=\frac{-3u(x_1) + 4u(x_2) -u(x_3)}{2h}$
.
Queste formule si ottengono sviluppando con Taylor e cercando un combinazione lineare opportuna in
questo modoIl miglior riferimento in italiano che ho sull'argomento è
Quarteroni-Matematica Numerica, dove c'è un intero capitolo sull'argomento.
Ora, visto che la condizione è sul bordo $x=0$, devi scrivere la corrispondente discretizzazione della prima riga:
Al bordo $x=0$ si ha che
$\frac{u_2 - u_0}{2h}= 1 $
da cui $u_0 = u_2 - 2h$.
La prima riga è perciò
$u_1 + \frac{u_2-u_0}{2h} + \frac{u_0 - 2 u_1 + u_2}{h^2}=\cos(x_1)$
Pertanto, ricordando il valore per $u_0$ appena ricavato, e che $\cos(x_1)=1$, si ha
$u_1 + \frac{2}{h^2} u_2 - \frac{2}{h} - 2 \frac{u_1}{h^2}=0$
Quindi la prima riga, finalmente, è
$u_1 (1 - \frac{2}{h^2}) + \frac{2}{h^2} u_2 = \frac{2}{h}$
Nel tuo caso, data la matrice del problema, cioè $D =A + B + I$, si ha che va imposto che
la prima riga della discretizzazione sia uguale a quanto ricavato sopra.
Questo può essere fatto in MatLab/Octave con il comando
- Codice:
D(1,1:2)=[1 - (2/h^2),2/h^2]
, che cambia le prime due entrate della matrice, e modificando il termine noto $b=\cos(x)$, con il comando
- Codice:
b(1)=2/h
.
Note: 1. il mio $\text{b}$ corrisponde al tuo $\text{rhs}$
2. Questo è un caso speciale perché il problema è lineare. Ma se l'equazione da azzerare fosse stata non-lineare, come poteva essere, ad esempio, $u''(x) + u(x) = cos(u(x))$, allora la $F(\vec{u})$ che risulta è un funzionale non lineare, e l'equazione $F(\vec{u})=0$ va risolta con un metodo di Newton/ Newton modificato.
Condizione di DirichletQui è molto più facile, va imposto che nell'ultima riga valga $u_m \approx u (x_m)=1$.
A tal fine, un modo (poco carino), è quello di azzerare l'ultima riga delle matrici di discretizzazione $A, B$ con il seguente comando:
- Codice:
A(m,m-1:m)=[0,0]; B(m,m-1:m)=[0,0]
Ora l'ultima riga si legge:
$u_m =cos(x_m)$
Poiché $\cos(x_m)$ è nullo, va cambiato con
- Codice:
b(m)=1;
Ora l'ultima riga è sistemata.
RisoluzioneIl problema è lineare, ed è, precisamente,
$D \vec{u}= \vec{b}$
e si risolve con uno dei metodi per la risoluzione di sistemi lineari (oppure col backslash di MatLab).
ImplementazioneUna possibile implementazione, facendo riferimento al tuo codice, potrebbe essere:
- Codice:
clc
clear all
close all
mrif=2^10+1;
m=mrif;
x=linspace(0,pi/2,m)';
h=(pi/2)/ (m-1);
%discretizzazione derivate
A = toeplitz(sparse([1,1],[1,2],[-2,1]/(h^2),1,m));
B = toeplitz(sparse(1,2,-1/(2*h),1,m), sparse(1,2,1/(2*h),1,m));
I=speye(m);
%% CONDIZIONI AL BORDO
%DIRICHLET al nodo a destra.
A(m,m-1:m)=[0,0];
B(m,m-1)=0; % AZZERO PER FAR SOPRAVVIVERE SOLO LA MATRICE IDENTITÀ NELL'ULTIMA RIGA
D=A+B+I;
% NEUMANN, risolte con nodo fantasma
D(1,1:2) = [1-2/h^2,2/h^2];
b=cos(x); %termine noto
b(1)=2/h;
b(m)=1;
F=@(u) D*u - b; %funzione da azzerare
u=D\b;
figure
plot(x,u,'b',x,sin(x),'r-o')
axis([0, pi/2,0,1])
legend('sol numerica','sol analtica')
%% Costruzione soluzione di riferimento
Urif=u;
mrange=2.^(4:8)+1;
count=0;
for m=mrange
count=count+1;
x=linspace(0,pi/2,m)';
h=(pi/2)/ (m-1);
A = toeplitz(sparse([1,1],[1,2],[-2,1]/(h^2),1,m));
B = toeplitz(sparse(1,2,-1/(2*h),1,m), sparse(1,2,1/(2*h),1,m));
I=speye(m);
%CONDIZIONI AL BORDO
A(m,m-1:m)=[0,0];
B(m,m-1)=0; % AZZERO PER FAR SOPRAVVIVERE SOLO LA MATRICE IDENTITÀ NELL'ULTIMA RIGA
D=A+B+I;
D(1,1:2) = [1-2/h^2,2/h^2];
b=cos(x);
b(1)=2/h;
b(m)=1;
u=D\b;
err(count) = norm(u-Urif(1:(mrif-1)/(m-1):mrif),inf)
end
figure
loglog(mrange,err,'o',mrange,mrange.^-2*(err(end)/mrange(end)^-2))
che produce i grafici seguenti (dove nel secondo è riportato il logaritmo dell'errore in norma infinito).