posso chiederti come è stata calcolata la jacobiana? inoltre possiamo sempre scrivere l'equazione differenziale in forma matriciale? perchè per forme più complicate non credo si possa fare . . .
Il metodo delle differenze finite si basa proprio sul fatto di discretizzare l'equazione differenziale in modo tale da ottenere un sistema di equazioni algebriche da risolvere. Tale sistema di equazione è per l'appunto scrivibile in forma matriciale. La soluzione di tale sistema è il vettore delle approssimazioni della soluzione valutata nei punti di discretizzazione.
Per quanto riguarda il calcolo della jacobiana, è una cosa abbastanza magica e a volte incasinatrice. le prime volte è possibile provare a derivare manualmente elemento per elemento per vedere con i propri occhi come funziona... dopodiché ci si può abbastanza fidare del fatto che con qualche accorgimento tutto funziona similmente al caso scalare. Bisogna prestare molta attenzione alle dimensioni degli oggetti e a eventuali trasposizioni di vettori. Innanzitutto se F è un vettore, JF sarà una matrice. Ad esempio: A*y ---> derivo rispetto a y ---> A y.^2 ---> derivo rispetto a y ---> diag(2*y) [è una matrice!]
clear close clc %metodo risoluzioni finite: guarda one note m=51; x=linspace(0,1,m); h=x(end)/(m-1); e=ones(m,1); %matrice che approssima derivata prima di y D1=spdiags([-e/(2*h),zeros(m,1),e/(2*h)],-1:1,m,m); %matrice che approssima derivata seconda di y D2=spdiags([e/(h^2),-2*e/(h^2),e/(h^2)],-1:1,m,m); %termine noto b=exp(2*x); %Condizioni di Dirichlet: y(0)=1 D1(1,:)=0; D2(1,:)=0; b(1)=1; %Condizioni di Dirichlet: y(1)=3 D1(m,:)=0; D2(m,:)=0; b(m)=9; %sistema di equazioni non lineari: Metodo di Newton F=@(y) D2*y-3*D1*y+y.^2-b; JF = @(y) D2 - 3*D1 + diag(2*y); %Jacobiana y0 = ones(m, 1) * 1.5; %primo tentativo: ho scelto 2 perchè è tra 1 e 3 tol = h^2 / 100; %tolleranza: o.d.g diviso per tanti 0 quanto l'ordine dell'equazione y = y0; %inizializzazione delta = -JF(y) \ F(y); % Newton-Raphson iteration while (norm(delta, inf) > tol) y = y + delta; delta = -JF(y) \ F(y); end y=y+delta; yrif=y; plot(x, yrif, '*') xlabel('x') ylabel('y')
@@diegorampi4711Forse è perché il tuo b è un vettore riga anziché un vettore colonna e Matlab fa casino! Poi tutto funzionerà :) Grazie mille per la domanda.
@@stefanospurio7641 Mi dispiace ma io non la conosco ;) come d'altronde nel 99% dei casi. Puoi provare ad utilizzare un risolutore di eq. diff. online oppure ad utilizzare qualche metodo di matlab (es: bvp4c)
sei il king
posso chiederti come è stata calcolata la jacobiana? inoltre possiamo sempre scrivere l'equazione differenziale in forma matriciale? perchè per forme più complicate non credo si possa fare . . .
Il metodo delle differenze finite si basa proprio sul fatto di discretizzare l'equazione differenziale in modo tale da ottenere un sistema di equazioni algebriche da risolvere.
Tale sistema di equazione è per l'appunto scrivibile in forma matriciale.
La soluzione di tale sistema è il vettore delle approssimazioni della soluzione valutata nei punti di discretizzazione.
Per quanto riguarda il calcolo della jacobiana, è una cosa abbastanza magica e a volte incasinatrice.
le prime volte è possibile provare a derivare manualmente elemento per elemento per vedere con i propri occhi come funziona...
dopodiché ci si può abbastanza fidare del fatto che con qualche accorgimento tutto funziona similmente al caso scalare.
Bisogna prestare molta attenzione alle dimensioni degli oggetti e a eventuali trasposizioni di vettori.
Innanzitutto se F è un vettore, JF sarà una matrice.
Ad esempio:
A*y ---> derivo rispetto a y ---> A
y.^2 ---> derivo rispetto a y ---> diag(2*y) [è una matrice!]
Buongiorno, il codice non gira (Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 2.122403e-19. )
clear
close
clc
%metodo risoluzioni finite: guarda one note
m=51;
x=linspace(0,1,m);
h=x(end)/(m-1);
e=ones(m,1);
%matrice che approssima derivata prima di y
D1=spdiags([-e/(2*h),zeros(m,1),e/(2*h)],-1:1,m,m);
%matrice che approssima derivata seconda di y
D2=spdiags([e/(h^2),-2*e/(h^2),e/(h^2)],-1:1,m,m);
%termine noto
b=exp(2*x);
%Condizioni di Dirichlet: y(0)=1
D1(1,:)=0;
D2(1,:)=0;
b(1)=1;
%Condizioni di Dirichlet: y(1)=3
D1(m,:)=0;
D2(m,:)=0;
b(m)=9;
%sistema di equazioni non lineari: Metodo di Newton
F=@(y) D2*y-3*D1*y+y.^2-b;
JF = @(y) D2 - 3*D1 + diag(2*y); %Jacobiana
y0 = ones(m, 1) * 1.5; %primo tentativo: ho scelto 2 perchè è tra 1 e 3
tol = h^2 / 100; %tolleranza: o.d.g diviso per tanti 0 quanto l'ordine dell'equazione
y = y0; %inizializzazione
delta = -JF(y) \ F(y);
% Newton-Raphson iteration
while (norm(delta, inf) > tol)
y = y + delta;
delta = -JF(y) \ F(y);
end
y=y+delta;
yrif=y;
plot(x, yrif, '*')
xlabel('x')
ylabel('y')
@@diegorampi4711Forse è perché il tuo b è un vettore riga anziché un vettore colonna e Matlab fa casino! Poi tutto funzionerà :) Grazie mille per la domanda.
@@math4every1_ confermo che scrivendo b come vettore colonna il codice gira e bene
@@math4every1_ si conosce la soluzione esatta per fare un confronto?
@@stefanospurio7641 Mi dispiace ma io non la conosco ;) come d'altronde nel 99% dei casi.
Puoi provare ad utilizzare un risolutore di eq. diff. online oppure ad utilizzare qualche metodo di matlab (es: bvp4c)