Differenze finite di secondo ordine

[stripdance4ever]

Ciao a tutti. Avrei un dubbio riguardante un codice in matlab per risolvere equazioni differenziali mediante differenze finite di secondo ordine; la consegna dice:
Si risolva il problema ai limiti | u''(x)+u'(x)+u(x)-cos(x) = 0 0 | u'(0) = 1
| u(pi/2) = 1
usando il metodo delle differenze finite del secondo ordine. Si mostri inoltre l'ordine del metodo mediante un grafico log-log dell'errore in norma infinito rispetto ad una soluzione di riferimento.

Questo è il codice:

mrif = 901;
m = mrif;
x = linspace(0,pi/2,m)';
h = (pi/2)/(m-1);
A = toeplitz(sparse([1,2],[1,1],[-2 1]/h^2,m,1));
B = toeplitz(sparse(2,1,-1/(2*h),m,1),sparse(2,1,1/(2*h),m,1));
I = speye(m);
A(1,2) = (2/h^2)+2; ???
A(m,m-1:m) = [0,0]; ???
B(1,2) = 0; ????
B(m,m-1) = 0; ????
D = A+B+I;
rhs = cos(x);
rhs(1) = rhs(1)-1+2/h ?????
rhs(end) = 1;
u = D\rhs;
urif = u;

mrange = [21 31 51 61 91 101 151];
counter = 1;
for m = mrange
x = linspace(0,pi/2,m)';
h = (pi/2)/(m-1);
A = toeplitz(sparse([1,2],[1,1],[-2 1]/h^2,m,1));
B = toeplitz(sparse(2,1,-1/(2*h),m,1),sparse(2,1,1/(2*h),m,1));
I = speye(m);
A(1,2) = 2/h^2; ????
A(m,m-1:m) = [0 0];
B(1,2) = 0;
B(m,m-1) = 0;
D = A+B+I;
rhs = cos(x);
rhs(1) = rhs(1)-1+2/h;
rhs(end) = 1;
u = D\rhs;
err(counter) = norm(u-urif(1:(mrif-1)/(m-1):end),Inf);
counter = counter+1;
end

loglog(mrange,err,'o',mrange,mrange.^-2*(err(end)/mrange(end)^-2))


Non capisco le parti in cui ci sono i punti di domanda, che sostanzialmente riguardano le condizioni al bordo. Se qualcuno riuscisse a darmi una mano. Grazie.

[feddy]

Ciao, benvenuto/a sul forum. Consiglio, come da [regolamento]regolamento[/regolamento], di scrivere in LaTeX. Vedi come si scrivono le formule: [formule][/formule]

il 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.


[size=150]Condizione di Neumann[/size]

Per 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

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 modo

Il 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

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

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.



[size=150]
Condizione di Dirichlet[/size]

Qui è 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:

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

b(m)=1;

Ora l'ultima riga è sistemata.


[size=150]Risoluzione[/size]

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

[size=150]Implementazione[/size]

Una possibile implementazione, facendo riferimento al tuo codice, potrebbe essere:

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

[stripdance4ever]

E' tutto chiaro tranne una cosa; non ho capito il passaggio per cui b=1.
Poi avrei un altro problema da proporti (adesso scriverò in LaTex ,scusa non lo sapevo) :

Si risolva il problema ai limiti :

$\{($(-del)/(delx)$)(($1+x$)*($(del)/(delx)$)*($u(x)$))=1,($u(0)=0$),($u(1)=0$):}$ con ${0
Si mostri inoltre l'ordine del metodo mediante un grafico loglog
dell'errore in norma infinito rispetto alla soluzione esatta
$(1+x)*u''(x)$+$u'(x)$=1

Spero di aver scritto correttamente in LaTex

[stripdance4ever]

Forse così è meglio:
$\{(-(del)/(delx)*((1+x)* (del)/(delx)*(u(x)))=1, (u(0)=0), (U(1)=0):}$

[feddy]

"s-o-p-h-i":non ho capito il passaggio per cui b=1.

Ho imposto che l'ultima riga si leggesse $u_m=1$. E' uguale a quello che hai fatto te, solo che è più chiaro problabilmente.

"s-o-p-h-i":Poi avrei un altro problema da proporti

Dovresti aprire un nuovo post a dire il vero, però visto che è molto simile fa lo stesso

Per imparare come scrivere giusto il testo, fai click destro sulla formula e premi su 'Show math as ' --> 'TeX commands'

Il problema è:

\( \begin{cases} -\frac{d}{dx}((1+x) \frac{d}{dx}u(x))=1 \\ u(0)=0 \\ u(1)=0 \end{cases} \)

Basta usare la derivata del prodotto e l'equazione differenziale da risolvere diventa semplicemente

\( (1+x) u''(x) + u'(x) -1 =0 \)

che è un semplice problema lineare (nota che anche se avessimo avuto $e^x$ invece che $(1+x)$ sarebbe stato lineare).
Indicando con $A$ e $B$ le matrici che discretizzano tramite differenze finite gli operatori differenziali si ha che il problema discretizzato si scrive come

$( (1+x) A+ B) \cdot u = \vec{1}$

[size=150]Condizioni al bordo[/size]

Senza azzerare componenti delle matrici, è sufficiente porre $u(1)=0$ e $u(m)=0$ nella funzione da azzerare.

Altrimenti, azzerando, si può fare così:
\( \bullet \) chiamiamo il termine noto

b=ones(m,1);

e imponiamo che

b(1)=0;b(m)=0;

. Tra un momento sarà chiaro perchè.

\( \bullet \) Imponiamo che la prima riga sia $u(1)=0$. Azzeriamo la prima riga di $B$ e nella prima riga di $A$ lasciamo solo che sopravviva l'entrata (1,1), cioè

A(1,1:2)=[1,0]; B(1,2)=0;

\( \bullet \) Analogamente a quanto fatto sopra, si azzera l'ultima riga di $B$ e si lascia che sopravviva solo l'entrata di $A$ in posizione (m,m). Così si ha che $(1+x_m)*u_m = b(m)=0$, da cui $u_m=0$.


[size=150]Risoluzione[/size]

Si può scrivere l'equazione raccogliendo

M=-diag(1+x)*D2-D1;

La funzione da azzerare è

$F(u) = M*\vec{u} -\vec{b};$

e il problema si risolve con

u=M\b

[size=150]Implementazione[/size]

Se usi MatLab devi sostituire il comando

count ++ 

con il comando

count = count +1

clear all
close all


mrange=[10,20,100,200];
count=0;

for m=mrange+1
    count++;
    h=1/(m-1);
    x=linspace(0,1,m)';
    A = toeplitz(sparse([1,2],[1,1],[-2,1]/(h^2),m,1)); 
    B = toeplitz(sparse(2,1,-1/(2*h),m,1), sparse(2,1,1/(2*h),m,1));
   
    %boundary cond. 
    A(1,1:2)=[1,0];
    A(m,m-1:m)=[0,1];
    B(1,1:2)=[0,0];
    B(m,m-1:m)=[0,0];
    
    b=ones(m,1);
    b(1)=0;
    b(m)=b(1);
    
    M=-diag(1+x)*A-B; %matrice del problema    
    u=M\b; 
    
    s=@(x) -x + (log(1 + x))/(log(2)); %sol analitica
    err(count)=norm(u-s(x),inf);
end

figure(1)
plot(x,s(x),'k',x,u,'o')
axis([0,1,0,0.09])
legend('sol analitica','sol numerica')
figure(2)
loglog(mrange,err,'*',mrange,(err(end))*(mrange/mrange(end)).^-2)
legend('errore','ordine 2')

che riproduce con il corretto ordine la soluzione analitica

[stripdance4ever]

Non ho capito come hai fatto a trovare la soluzione analitica,che formula ci sta dietro?
comunque grazie,sei molto bravo,sicuramente ti proporrò altri problemi
probabilmente hai fatto o stai facendo il mio stesso corso di studio.

[feddy]

"s-o-p-h-i":Non ho capito come hai fatto a trovare la soluzione analitica,che formula ci sta dietro?

Nessuna formula, da

$\frac{d}{dx} ((1+x) u'(x)) = -1$

integri successivamente due volte e con le condizioni date ricavi la costante.

In pratica, integrando una volta trovi $[u'(x)= \frac{C- x}{1+x}]$.

Integrando una seconda:

$[u(x)=(C+1) \log(x+1) - x]$

Imponendo le condizioni (che non possono essere contraddittorie) si vede subito che

$C+1=\frac{1}{\log(2)}$

.

In definitiva la soluzione è

$u(x)=\frac{1}{\log(2)} \log(x+1) - x$