Equazione di poisson

[prezzemolina86]

Questo è il testo dell'esercizio...

Utilizzando l’approssimazione alle differenze finite con la regola dei 5 punti,
risolvere il problema al contorno:
Δu = -4(pi^2)u;
nel dominio quadrato [-1; 1]^2 e condizioni al contorno nulle. Si utilizzi il
metodo iterativo di Jacobi per la risoluzione del sistema lineare associato.
Si faccia il confronto con la soluzione esatta e si determini numericamente
l’ordine di approssimazione. Si stabilisca inoltre la convergenza del metodo
iterativo considerato. Si conservino i dati ottenuti in files :f ig di MATLAB.

giù c'è la function[size=150].il mio problema è che non capito come si passa la u[/size]...non
è quello che cerco??help!thanks!

function [x,y,u,uexact,errore,nit,dx,B,D,U0,L0]=poisson_5punti_Jacobi(L1,L2,M,BCd,BCu,BCr,BCl,f,toll)

x=linspace(L1,L2,M+2);
y=linspace(L1,L2,M+2);

dx=(L2-L1)/(M+1);

nmax=2000;

Uold=zeros(M^2,1);
U=zeros(M^2,1);
b=zeros(1,M^2);
err=zeros(1,nmax);

C=toeplitz([-4,1,zeros(1,M-2)]);
B1=blkdiag(C);

for k=1:M-1
B1=blkdiag(B1,C);
end

B2=toeplitz([zeros(1,M),1,zeros(1,M^2-(M+1))]);
B=B1+B2;

for i=1:M^2
D(i,i)=B(i,i);
end

for i=2:M^2
for j=1:M^2
if j L0(i,j)=-B(i,j);
end
end
end
for i=1:M^2
for j=1:M^2
if j>i
U0(i,j)=-B(i,j);
end
end
end
U0=[U0;zeros(1,M^2)];
L0=[L0,zeros(1,M^2)'];

P=D;
N=L0+U0;
H=inv(P)*N;

u(1,:)=feval(BCl,y);
u(M+2,:)=feval(BCr,y);
u(:,1)=feval(BCd,x);
u(:,M+2)=feval(BCu,x);

for i=1:M
for j=1:M
b((j-1)*M+i)=dx^2*feval(f,x(i+1),y(j+1));
end
end

b(1:M:end)=b(1:M:end)-u(1,2:end-1);
b(M:M:end)=b(M:M:end)-u(end,2:end-1);

b(1:M)=b(1:M)-u(2:end-1,1)';
b(end-M+1:end)=b(end-M+1:end)-u(2:end-1,end)';

v=inv(P)*b';

for n=1:nmax
U=H*Uold+v;
err(n)=max(abs(U-Uold));
nit=n;
if err(n) break;
else
Uold=U;
end
end

for i=1:M
for j=1:M
u(i+1,j+1)=U((j-1)*M+i);
end
end

%Equazione di Laplace
% for i=1:M+2
% for j=1:M+2
% uexact(i,j)=y(j)/((1+x(i))^2+y(j)^2);
% end
% end

%Equazione del compito
for i=1:M+2
for j=1:M+2
uexact(i,j)=sin(sqrt(2)*pi*(x(i)+y(j)));
end
end
%Equazione con f=x^2+y^2
%for i=1:M+2
% for j=1:M+2
% uexact(i,j)=exp(pi*x(i))*sin(pi*y(j))+(1/2)*((x(i)*y(j))^2);
% end
%end

errore=max(max(abs(u-uexact)));

[hamming_burst]

Ciao,
prova a mettere il codice nei tag così si mantiene l'identazione e si capisce qualcosa, scusa ma così è inguardabile

[prezzemolina86]

testo esercizio...

Utilizzando l’approssimazione alle differenze finite con la regola dei 5 punti,
risolvere il problema al contorno:
Δu = -4(pi^2)u;
nel dominio quadrato [-1; 1]2 e condizioni al contorno nulle. Si utilizzi il
metodo iterativo di Jacobi per la risoluzione del sistema lineare associato.
Si faccia il confronto con la soluzione esatta e si determini numericamente
l’ordine di approssimazione. Si stabilisca inoltre la convergenza del metodo
iterativo considerato. Si conservino i dati ottenuti in files :f ig di MATLAB.

spero che adesso si vedrà meglio...ricordo il problema...[size=150]come passo la u?[/size]pensavo che qst function
servisse ha trovare la u!

function [x,y,u,uexact,errore,nit,dx,B,D,U0,L0]=poisson_5punti_Jacobi(L1,L2,M,BCd,BCu,BCr,BCl,f,toll)

x=linspace(L1,L2,M+2);
y=linspace(L1,L2,M+2);

dx=(L2-L1)/(M+1);

nmax=2000;

Uold=zeros(M^2,1);
U=zeros(M^2,1);
b=zeros(1,M^2);
err=zeros(1,nmax);

C=toeplitz([-4,1,zeros(1,M-2)]);
B1=blkdiag(C);

for k=1:M-1
    B1=blkdiag(B1,C);
end

B2=toeplitz([zeros(1,M),1,zeros(1,M^2-(M+1))]);
B=B1+B2;

for i=1:M^2
    D(i,i)=B(i,i);
end

for i=2:M^2
    for j=1:M^2
        if j<i
            L0(i,j)=-B(i,j);
        end
    end
end
for i=1:M^2
    for j=1:M^2
        if j>i
            U0(i,j)=-B(i,j);
        end
    end
end
U0=[U0;zeros(1,M^2)];
L0=[L0,zeros(1,M^2)'];

P=D;
N=L0+U0;
H=inv(P)*N;

u(1,:)=feval(BCl,y);
u(M+2,:)=feval(BCr,y);
u(:,1)=feval(BCd,x);
u(:,M+2)=feval(BCu,x);

for i=1:M
    for j=1:M
        b((j-1)*M+i)=dx^2*feval(f,x(i+1),y(j+1));
    end
end

b(1:M:end)=b(1:M:end)-u(1,2:end-1);
b(M:M:end)=b(M:M:end)-u(end,2:end-1);

b(1:M)=b(1:M)-u(2:end-1,1)';
b(end-M+1:end)=b(end-M+1:end)-u(2:end-1,end)';

v=inv(P)*b';

for n=1:nmax
    U=H*Uold+v;
    err(n)=max(abs(U-Uold));
    nit=n;
    if err(n)<toll
        break;
    else
        Uold=U;
    end
end

for i=1:M
    for j=1:M
        u(i+1,j+1)=U((j-1)*M+i);
    end
end

%Equazione di Laplace
% for i=1:M+2
%     for j=1:M+2
%         uexact(i,j)=y(j)/((1+x(i))^2+y(j)^2);
%     end
% end

%Equazione del compito
 for i=1:M+2
     for j=1:M+2
         uexact(i,j)=sin(sqrt(2)*pi*(x(i)+y(j)));
     end
 end
%Equazione con f=x^2+y^2
%for i=1:M+2
%   for j=1:M+2
%       uexact(i,j)=exp(pi*x(i))*sin(pi*y(j))+(1/2)*((x(i)*y(j))^2);
%   end
%end

errore=max(max(abs(u-uexact)));

[lukul]

Ciao,

prima di tutto, data la funzione che hai definito nella prima riga

function [x,y,u,uexact,errore,nit,dx,B,D,U0,L0]=poisson_5punti_Jacobi(L1,L2,M,BCd,BCu,BCr,BCl,f,toll)

cerca di spiegare ciascuna variabile di input e di output a cosa corrisponde (immagino che toll indichi la tolleranza, ma il resto lo ignoro).
Il resto a seguire...

[prezzemolina86]

testo esercizio...

Utilizzando l’approssimazione alle differenze finite con la regola dei 5 punti,
risolvere il problema al contorno:
Δu = -4(pi^2)u;
nel dominio quadrato [-1; 1]2 e condizioni al contorno nulle. Si utilizzi il
metodo iterativo di Jacobi per la risoluzione del sistema lineare associato.
Si faccia il confronto con la soluzione esatta e si determini numericamente
l’ordine di approssimazione. Si stabilisca inoltre la convergenza del metodo
iterativo considerato. Si conservino i dati ottenuti in files :f ig di matlab.

spero che adesso si vedrà meglio...ricordo il problema...come passo la u?pensavo che qst function
servisse ha trovare la u!

Codice:
function [x,y,u,uexact,errore,nit,dx,B,D,U0,L0]=poisson_5punti_Jacobi(L1,L2,M,BCd,BCu,BCr,BCl,f,toll)
%L1,L2 sono gli estremi del mio dominio
%M è il numero di punti in cui suddivido il dominio
%BCd,BCu,BCr,BCl sono rispettivamente condizioni al contorno in basso, in alto, a destra ed a sinistra
%f sarebbe la f!cioè -4*(pi^2)*u   si tratta proprio di questa u!

%creo la griglia
x=linspace(L1,L2,M+2);
y=linspace(L1,L2,M+2);

dx=(L2-L1)/(M+1);

nmax=2000;

%creo prima vettori e matrici per poi cosruire il sistema da risolvere
Uold=zeros(M^2,1);
U=zeros(M^2,1);
b=zeros(1,M^2);
err=zeros(1,nmax);

C=toeplitz([-4,1,zeros(1,M-2)]);
B1=blkdiag(C);

for k=1:M-1
    B1=blkdiag(B1,C);
end

B2=toeplitz([zeros(1,M),1,zeros(1,M^2-(M+1))]);
B=B1+B2;

%poichè voglio usare il metodo iterativo di jacobi per risolvere
%il sistema calcolo la matrice diaginale, triangolare inferiore e superiore
%la soluzione sarà x^k+1=D^-1*(L0+U0)x^k+(D^-1)b
for i=1:M^2
    D(i,i)=B(i,i);
end

for i=2:M^2
    for j=1:M^2
        if j<i
            L0(i,j)=-B(i,j);
        end
    end
end
for i=1:M^2
    for j=1:M^2
        if j>i
            U0(i,j)=-B(i,j);
        end
    end
end
U0=[U0;zeros(1,M^2)];
L0=[L0,zeros(1,M^2)'];

P=D;
N=L0+U0;
H=inv(P)*N;

u(1,:)=feval(BCl,y);
u(M+2,:)=feval(BCr,y);
u(:,1)=feval(BCd,x);
u(:,M+2)=feval(BCu,x);

for i=1:M
    for j=1:M
        b((j-1)*M+i)=dx^2*feval(f,x(i+1),y(j+1));
    end
end

b(1:M:end)=b(1:M:end)-u(1,2:end-1);
b(M:M:end)=b(M:M:end)-u(end,2:end-1);

b(1:M)=b(1:M)-u(2:end-1,1)';
b(end-M+1:end)=b(end-M+1:end)-u(2:end-1,end)';

v=inv(P)*b';

for n=1:nmax
    U=H*Uold+v;
    err(n)=max(abs(U-Uold));
    nit=n;
    if err(n)<toll
        break;
    else
        Uold=U;
    end
end

for i=1:M
    for j=1:M
        u(i+1,j+1)=U((j-1)*M+i);
    end
end


%cosruisco la soluzione esatta per poi fare il confronto con la 
%soluzione numerica
 for i=1:M+2
     for j=1:M+2
         uexact(i,j)=sin(sqrt(2)*pi*(x(i)+y(j)));
     end
 end


errore=max(max(abs(u-uexact)));

[lukul]

L’ideale per un utente esterno ma anche per lo stesso programmatore che diviene, in un secondo momento, utente di sé stesso, èquello di poter disporre di una struttura programmativa cosiffatta:

Routine
dedicata alla risoluzione del problema numerico specifico (ad esempio azzeramento di funzione tramite metodo dicotomico) in grado di ricevere via lista il nome della funzione da azzerare, l’intervallo di incertezza iniziale [a,b], il numero massimo di iterazioni da effettuare MAX_ITER, la massima ampiezza finale accettabile dell’intervallo di incertezza δ.
Tale routine dovrebbe restituire il punto, c, candidato alla soluzione in accordo con la struttura dell’algoritmo numerico.
L’utente, potendo disporre della routine summenzionata non avrebbe altro che da fornire il nome della funzione da azzerare e gli estremi dell’intervallo di incertezza a, bsoddisfacenti la condizione di Bolzano: f(a)·f(b) ≤0nonchécome indicato in precedenza i due parametri MAX_ITERe δ.

L'operazione di passare ad una routine il nome di un'altra routine (nel nostro caso funzione) è definito nel linguaggio tecnico "passaqggio di un puntatore di funzione".

Per passare via lista ad una routine di Matlab™il nome di una funzione occorre anteporre al nome della funzione stessa il simbolo: @.
Ad esempio se la funzione si chiama MyFunc, il suo puntatore viene passato alla routine WorkWithFunc tramite l’istruzione:
ret= WorkWithFunc(@MyFunc, ... )
-----
Supponiamo che la funzione MyFunc sia cosìstrutturata:
functionf = MyFunc(x)
f = x ^ 2 * sin(x) –cos(3. * x) + 1.;
-----

Nel tuo caso dovresti creare una function per.es. MyFunc contenente la funzione di interesse (la tua u) e poi effettuare la chiamata dalla function che hai su riportato mediante la sintassi su vista ( nome_variabile = Nome_function_da_cui_effettui_la_chiamata(@nome_function_che_vuoi_passare)

[lukul]

Cmque, onde evitare confusione, ti consiglio di cercare su internet qualcosa sui "puntatori a funzione" usati nel matlab. Tra l'altro sono gli unici puntatori di cui il matlab fà uso. Mi auguro di essere stato di aiuto. A disposizione per ogni altra domanda. Ciao...

[prezzemolina86]

grazie,cmq non era qst il punto. nel senso ke @(x) o @(y) li uso già quando passo le condizioni al contorno...mi sono spiegata male. non so COSA passare per u e non COME...
cmq ho provato a fare l esercizio con un problema di laplace per togliermi momentaneamente il problema della u.precisamente
$ Deltau=0 $
$ u(x,0)=0 (BCdown) $ $ u(x,1)=1/((1+x^2)+1) (BCup) $
$ u(0,y)=y/(1+y^2) (BCleft) $ $ u(1,y)=y/(4+y^2) (BCright) $
nel dominio [0,1]^2;

Ricordo il testo:

Utilizzando l’approssimazione alle differenze finite con la regola dei 5 punti,
risolvere il problema al contorno:
Δu = -4(pi^2)u;
nel dominio quadrato [-1; 1]^2 e condizioni al contorno nulle. Si utilizzi il
metodo iterativo di Jacobi per la risoluzione del sistema lineare associato.
Si faccia il confronto con la soluzione esatta e si determini numericamente
l’ordine di approssimazione. Si stabilisca inoltre la convergenza del metodo
iterativo considerato. Si conservino i dati ottenuti in files :f ig di matlab.

ho fatto il confronto tra la soluzione esatta e quella numerica.coincidono.quindi ok
Ora:si determini l ordine di approssimazione.io ero convinta ke l ordine di approssimazione
si facesse vedere mettendo su un grafico sull asse delle x il deltax e sull'asse delle ascisse
l'errore.ma penso che sia sbagliato.
negli appunti ho scritto

La consistenza è l'ordine di approssimazione.diremo che l'ordine di approssimazione è p se
$ u^(n+1)_(Deltax)=A_(Deltax)u^n_(Deltax)+k^n_(Deltax)+O(Deltax)^(p+r) $ dove questa u
è la soluzione esatta calcolata nei punti della griglia.

qualcuno mi sa spiegare quello che dovrei fare?GRAZIE MILLE.se non capisco il testo
non ho dove andare!