Vinnaren i pepparkakshustävlingen!
  1. Vetenskap & humaniora
  2. Fysik, matematik och teknologi
  3. Matematiska och naturvetenskapliga uppgifter

Lösning av PDE i Matlab, "diskretisering" (FDM)

Svara
2013-11-11, 01:50
  #1
esunix
Medlem
Hej,

Givet är ett "BVP", med

Laplace T = 0

med Dirichlet villkor
T(0,y) = 300, T(4,y) = 600 0<y<2

och Neumanvillkor
dT/dx(x,0)=0, dT/dy=(x,2)=0 0<x<4

I problemet skall man först "diskretisera" (med finita differensmetoden) området till en kvadratisk yta med konstant steglängd h=0.2 och i både x- och y-led samt med N=19, M=11. Hur är kvadratisk möjligt egentligen då vi har en rektaningel?

Vidare, på vilket sätt kan jag anpassa följande matlabkod till problemet ovan? Jag har lite problem att förstå hur många okända ekvationer ovanstående problem egentligen ger.
Matlabkoden löser ett laplaceekvationen i ett kvadratisk område med sidan 1 och boundary condition: u=0


Kod: 
clear all;

N = 100;

hx=1/(N+1); hy=1/(N+1);

x=0:hx:1; y=0:hy:1;
d=4*ones(N*N,1);
d_1=-ones(N*N,1);
d1=-ones(N*N,1);
for i=N:N:N*N-1
    d_1(i)=0;
    d1(i+1)=0;
end

d_2=-ones(N*N,1);
d2=-ones(N*N,1);
a = spdiags([d_2 d_1 d d1 d2], [-N, -1:1,N],N*N,N*N);

f=ones(N*N,1);
uinner=a\f;
Uinner=zeros(N,N);

for i=1:N
    Uinner(i,1:N)=uinner((i-1)*N+[1:N]);
end
U=zeros(N+2,N+2);
for i = 1:N
    U(i+1,2:N+1)=Uinner(i,1:N);
end
mesh(x,y,U)
title('Solution of \Laplace u = 1 on the unit square')
    xlabel('x')
    ylabel('y')
    zlabel('u(x,y)')

Tusen tack på förhand!
Citera
2013-11-11, 11:01
  #2
trekantshatt
Medlem
De första raderna är ett dåligt kommenterat sätt att implementera FD schemat, sen vette fan vad det är för sätt randvillkoren hanteras (det ämnet glossades över som satan i kursen jag gick så vet inte hur man gör). Det kan hjälpa att kika här http://www.mathworks.se/help/matlab/...al-arrays.html .
__________________
Senast redigerad av trekantshatt 2013-11-11 kl. 11:04.
Citera
2013-11-11, 14:02
  #3
esunix
Medlem
Citat:
Ursprungligen postat av trekantshatt
De första raderna är ett dåligt kommenterat sätt att implementera FD schemat, sen vette fan vad det är för sätt randvillkoren hanteras (det ämnet glossades över som satan i kursen jag gick så vet inte hur man gör). Det kan hjälpa att kika här http://www.mathworks.se/help/matlab/...al-arrays.html .

den där 4ones(NN,1) tror jag är koeficienten då man Diskretiserar U i Position(i,j)
Tack för länken vet dock inte om det hjälper Direkt. Någon färsk numerisk Expert som Kan?
Citera
2013-11-11, 18:32
  #4
trekantshatt
Medlem
Orkar tyvärr inte tränga igenom linjära indexeringen men för dirichletvillkoren verkar de ju ha använt en slags penalty barrier metod kombinerat med forcerade värden? Antar att man kan implementera von neumannvillkor (med eller utan ghostpoints och interpolering) med en liknande approach med antingen ett skevt differansschema eller något annat längs ränderna.
Citera
2013-11-12, 00:34
  #5
esunix
Medlem
Tror du jag kan först "diskretisera" området på på ett sådant sätt att man har

(M-1) okända punkter i x-led och N-1 i y-led sådan att man får totalt (M-1)*(N-1) okända ekvationer i uppgiften ifråga?

Sedan ha för Dirichletvillkor: T_0,j=300 för j = 1,2,3...N och T_4,j = 600 för j = 1,2,3...N
För neumann: ?? ingen aning.
Ingen som vet?
Numeriska metoder är verkligen inte min grej... :S
Citera
2013-11-12, 09:25
  #6
-taki-
Medlem
-taki-s avatar
Jag orkar inte gå igenom din kod men jag brukar göra som på bilden vid dessa randvillkor.

Vid Dirichlet går jag fram med diskretiseringen till delta-x från randen med mittnoden (med diskretiseringen för andraderivatan som jag skrivit i bilden). Då hamnar en nod på randen men dess värde är ju känt så det värdet bakar jag in i differensekvationen.

Vid Neuman brukar jag dra diskretiseringen så att mittnoden hamnar på randen. Uttrycket under till höger beskriver det kända flödet och används sedan för att eliminera Tw. På samma sätt som ovan så bakas detta sedan in i differensekvationen som en känd kvantitet.

http://i1061.photobucket.com/albums/...ps28122282.jpg
__________________
Senast redigerad av -taki- 2013-11-12 kl. 09:27.
Citera
2013-11-12, 12:40
  #7
esunix
Medlem
Citat:
Ursprungligen postat av -taki-
Jag orkar inte gå igenom din kod men jag brukar göra som på bilden vid dessa randvillkor.

Vid Dirichlet går jag fram med diskretiseringen till delta-x från randen med mittnoden (med diskretiseringen för andraderivatan som jag skrivit i bilden). Då hamnar en nod på randen men dess värde är ju känt så det värdet bakar jag in i differensekvationen.

Vid Neuman brukar jag dra diskretiseringen så att mittnoden hamnar på randen. Uttrycket under till höger beskriver det kända flödet och används sedan för att eliminera Tw. På samma sätt som ovan så bakas detta sedan in i differensekvationen som en känd kvantitet.

http://i1061.photobucket.com/albums/...ps28122282.jpg

Skulle du kunna förklara mer om Neumann? Är T_e den punkten innanför randen då?
Citera
2013-11-12, 19:18
  #8
-taki-
Medlem
-taki-s avatar
Citat:
Ursprungligen postat av esunix
Skulle du kunna förklara mer om Neumann? Är T_e den punkten innanför randen då?

Ja, avståndet delta-x innanför. Mittpunkten i "molekylen", Tp, ligger på randen för dess värde vill man ju också åt.

Med ekvationen nederst till höger kan man sedan eliminera Tw.

Sedan finns ju den motsatta randen. Där hamnar Te utanför och ska elimineras.

Den molekyl jag skissat är för en endimensionell variant av diffekvationen i fråga. Vid 2D får man ju en Tn och Ts (Tnord och Tsyd) också i molekylen. d2T/dy2 = (Tn-2Tp+Ts)/dy2.

Du behöver alltså eliminera Tn för molekylerna där Tp ligger på "översta", norra, randen och Ts där Tp ligger på södra randen.

För alla noder där Tp ligger ett delta-x in från randen ser differensekvationerna identiska ut.
Åtta set av differensekvationer behöver specialskrivas. Det är de fyra i var sitt hörn och de som ligger med Tp på respektive rand.

Lite skrivjobb innan programmeringen blir det ...
Citera
2013-11-12, 22:46
  #9
esunix
Medlem
Nu har jag skapat "grids"

På x-axeln har jag ansatt i = 0, 1, 2,... M

y-axeln: j = 1, 2, 3,... N

dvs att på y-axeln startar man indexet från j = 1

Därmed av differensekvationen får jag:

4T_(i,j) + T_(i-1,j) + T_(i+1,j) + T_(i,j-1) + T_(i,j+1) = 0

Dirichletvillkor: T_(0,j) = 300 då j = 1,2,3...N och T_(M,j) = 600 då j = 1,2,3...N

För Neumann har jag skrivit:

(T_(i,0)-T_(i,2))/2h = 0 för i = 1,2,3....M-1

(T_(i,N+1)-T_(i,N-1))/2h = 0 för i = 1,2,3... M-1

Är detta rätt?
Undrar
Citera
2013-11-14, 16:06
  #10
-taki-
Medlem
-taki-s avatar
Citat:
Ursprungligen postat av esunix
4T_(i,j) + T_(i-1,j) + T_(i+1,j) + T_(i,j-1) + T_(i,j+1) = 0

Ska vara minustecken framför alla termer utom 4T:

4T_(i,j) - T_(i-1,j) - T_(i+1,j) - T_(i,j-1) - T_(i,j+1) = 0

i går väl lämpligast mellan 0 och M+1 och du sätter upp differensekvationen (mittnoden) mellan 1 <= i <= M. Noderna med M=0 och M+1 har kända värden. Då får du 19 noder med obekanta och h=0.2. Precis som uppgiften kräver.

Dirichletvillkor: T_(0,j) = 300 och T_(M+1,j) = 600 då j = 1,2,3...N


För Neumann då:

(T_(i,0)-T_(i,2))/2h = 0 och (T_(i,N+1)-T_(i,N-1))/2h = 0 för i = 1,2,3... M
Citera
Svara

Stöd Flashback

Flashback finansieras genom donationer från våra medlemmar och besökare. Det är med hjälp av dig vi kan fortsätta erbjuda en fri samhällsdebatt. Tack för ditt stöd!

Stöd Flashback