Differenze finite di secondo ordine

Messaggioda s-o-p-h-i » 03/06/2018, 21:29

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<x<pi/2
| 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. <Non capisco come si impostano.
Se qualcuno riuscisse a darmi una mano. Grazie.
s-o-p-h-i
Starting Member
Starting Member
 
Messaggio: 1 di 8
Iscritto il: 03/06/2018, 19:29

Re: Differenze finite di secondo ordine

Messaggioda feddy » 04/06/2018, 10:59

Ciao, benvenuto/a sul forum. Consiglio, come da Regolamento, di scrivere in LaTeX. Vedi come si scrivono le 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.


Condizione di Neumann

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
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 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
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 Dirichlet


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


Risoluzione

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

Implementazione

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

Immagine


Immagine
Avatar utente
feddy
Moderatore
Moderatore
 
Messaggio: 1888 di 5934
Iscritto il: 26/06/2016, 00:25
Località: SISSA

Re: Differenze finite di secondo ordine

Messaggioda s-o-p-h-i » 05/06/2018, 17:18

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<x<1}$

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 ](*,) :smt087
s-o-p-h-i
Starting Member
Starting Member
 
Messaggio: 2 di 8
Iscritto il: 03/06/2018, 19:29

Re: Differenze finite di secondo ordine

Messaggioda s-o-p-h-i » 05/06/2018, 17:33

Forse così è meglio:
$\{(-(del)/(delx)*((1+x)* (del)/(delx)*(u(x)))=1, (u(0)=0), (U(1)=0):}$
s-o-p-h-i
Starting Member
Starting Member
 
Messaggio: 3 di 8
Iscritto il: 03/06/2018, 19:29

Re: Differenze finite di secondo ordine

Messaggioda feddy » 05/06/2018, 18:04

s-o-p-h-i ha scritto: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 ha scritto:Poi avrei un altro problema da proporti


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

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}$


Condizioni al bordo

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
Codice:
b=ones(m,1);
e imponiamo che
Codice:
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è
Codice:
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$.


Risoluzione

Si può scrivere l'equazione raccogliendo
Codice:
M=-diag(1+x)*D2-D1;


La funzione da azzerare è
$F(u) = M*\vec{u} -\vec{b};$


e il problema si risolve con
Codice:
u=M\b



Implementazione

Se usi MatLab devi sostituire il comando
Codice:
count ++
con il comando
Codice:
count = count +1


Codice:
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

Immagine
Avatar utente
feddy
Moderatore
Moderatore
 
Messaggio: 1895 di 5934
Iscritto il: 26/06/2016, 00:25
Località: SISSA

Re: Differenze finite di secondo ordine

Messaggioda s-o-p-h-i » 06/06/2018, 13:24

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 :D
probabilmente hai fatto o stai facendo il mio stesso corso di studio.
s-o-p-h-i
Starting Member
Starting Member
 
Messaggio: 4 di 8
Iscritto il: 03/06/2018, 19:29

Re: Differenze finite di secondo ordine

Messaggioda feddy » 06/06/2018, 18:23

s-o-p-h-i ha scritto: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$
Avatar utente
feddy
Moderatore
Moderatore
 
Messaggio: 1896 di 5934
Iscritto il: 26/06/2016, 00:25
Località: SISSA


Torna a Analisi Numerica e Ricerca Operativa

Chi c’è in linea

Visitano il forum: Nessuno e 1 ospite