Problema con le condizioni di Neumann su Matlab

Salve a tutti,
ho un problema per l'implementazioni su Matlab di questo problema, ove la soluzione è circa 0.208, ma non ho idea di come fare le ho provate in tutti i modi ma non sono capace di risolvere Neumann con differenze centrate. Di sotto posto la traccia dell' esercizio sperando che qualcuno di voi riesca a darmi una mano.

Si consideri il problema:
u''(x) + sin(x) = 0 su un intervallo [0,2*pi]

con condizioni
u(0)=0
u'(2*pi) = 1

Si risolva con il metodo delle differenze finite, usando differenze centrate per l'approssimazione della condizione al bordo di Neumann.
Il valore della soluzione approssimata nel punto 2\pi ottenuta scegliendo una partizione uniforme dell'intervallo di integrazione in M=10 intervalli.

Al bordo $x=2 \pi $ devi usare un nodo fantasma, o "ghost node", in questo modo.
Supponi di discretizzare $[0, 2\pi]$ con $m$ nodi equispaziati. Allora, puoi discretizzare la condizione al bordo $u'(2 \pi )=1$ con $\frac{u_{m+1}-u_{m-1}}{2h}=1$. Da questa espressione ti ricavi (subito) quanto vale $u_{m+1}$. A questo punto, torni all'equazione discretizzata: $\frac{u_{i+1} - 2*u_i + u_{i-1}}{h^2} + \sin(x_i)=0$

Sostituisci qui dentro, per $i=m$, il valore di $u_{m+1}$ calcolato come descritto sopra: avrai così ottenuto l'ultima riga del tuo problema, e quindi dovrai modificare il sistema lineare corrispondente.

P.S.: Nota che con il ghost node si usa effettivamente un'approssimazione della derivata prima di ordine due, in accordo con quanto richiesto dal problema.

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.

Implementazione
Questo è 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!