$ { ( y''(x) + (y'(x))/x =cos(y(x))),( y'(0)=0 ),( y(0)=1 ):} , x \in (0,1]$
Lo risolvo con differenze finite. Poiché ho condizioni di Neumann sul primo nodo e sempre sullo stesso ho una di Dirichlet, ho pensato di porre $u_1=1$, in quanto $y(0)=1$, e poi utilizzare un ghost node $u_0$ per discretizzare la derivata prima.
In sostanza, detto $h$ il passo di discretizzazione, ho $y'(0)=(u_2-u_0) / (2h) = 0$, da cui $u_0=u_2$.
Sostituendo nell'equazione, correggetemi se sbaglio, trovo che la prima riga del sistema deve risultare uguale a
$2/h^2 u_2 - 2/h^2 -cos(1)=0$
.A tal fine, detta $D_1$ la matrice di discretizzazione della derivata prima, $D_2$ quella della seconda e $I$ l'identità, ho che il problema discretizzato equivale a risolvere l'equazione $F(u)= D_2 * u + 1/x * (D_1*u) - cos(u) =0$, dove in realtà, per non rendere pesante la notazione, qui $u=vecu$.
Evidentemente, c'è un problema con quel zero al denominatore: per aggirarlo ho considerato un linspace che parte da $0+eps$ e va fino a $1$.
Allego il codice del problema (in MatLab)
- Codice:
m=100;
h=(1-eps)/(m-1);
x=linspace(0+eps,1,m)';
%matrici di discretizzazione + boundary conditions
D2 = toeplitz(sparse([1,1],[1,2],[-2,1]/(h^2),1,m));
D1 = toeplitz(sparse(1,2,-1/(2*h),1,m), sparse(1,2,1/(2*h),1,m));
I=speye(m);
D2(1,1:2)=[0,2/(h^2)];
D1(1,2)=0;
I(1,1)=0;
b=zeros(m,1);
b(1,1)= 1 -2/h^2 - cos(1);
%funzione
Fun=@(u) D2*u + spdiags(1./x,0,m,m)*(D1*u) - cos(I*u) + b;
Jac=@(u) D2 + spdiags(1./x,0,m,m)*D1 + diag(sin(I*u));
%newton
u0=ones(m,1); %guess iniziale
tol=1e-10;
maxit=200;
u=u0;
delta=-Jac(u)\Fun(u);
i=0;
%Newton per sistemi NON-LINEARI
while(norm(delta,inf)>tol && i<maxit)
i=i+1;
u=u+delta;
delta=-Jac(u) \ Fun(u);
end
u=u+delta;
plot(x,u,'*')
Il plot che ne risulta è il seguente: dalle condizioni al contorno mi aspetterei tangente orizzontale, invece qui è proprio il contrario. Sbaglio a applicare differenze finite oppure c'è un errore nell'implementazione?
Grazie a chiunque mi possa illuminare