IMPLEMENTASI MATLAB UNTUK MENYELESAIKAN MASALAH SYARAT BATAS PERSAMAAN DIFFERENSIAL POISSON DAN LAPLACE 2D. La Ode Muhammad Umar RRR 1)

Paradigma, Vol. 14 No. 2 Agustus 2010 hlm. 151–170

IMPLEMENTASI MATLAB UNTUK MENYELESAIKAN MASALAH SYARAT BATAS PERSAMAAN DIFFERENSIAL POISSON DAN LAPLACE 2D La Ode Muhammad Umar RRR1) 1)

Jurusan Matematika, FMIPA, Universitas Haluoleo Kendari, 93231

ABSTRAK Pada tulisan ini, metode beda hingga dan metode iterasi Succesive Over Relaxation (SOR) digunakan untuk menyelesaikan persamaan differensial Poisson dan Laplace 2D pada domain segiempat. Metode ini tidak hanya menghasilkan akurasi numerik yang sangat baik tetapi juga sangat efisien. Lebih lanjut, diberikan program komputer dalam Matlab. Kata Kunci: Persamaan Poisson dan Laplace, metode beda hingga, SOR, Matlab

ABSTRACT In this paper, the finite difference methods and succesive over relaxation (SOR) are used to determine the solutions of Poisson and Laplace equation on a rectangular domain. These methods not only preserve the accuracy but also provide the efficiency. Moreover, computer program are presented in Matlab Keywords: Poisson and Laplace equations, finite difference, SOR, Matlab Diterima : 16 Juni 2010 Disetujui untuk dipublikasikan : 2 Agustus 2010

1. Pendahuluan Persamaan differensial Poisson dan Laplace 2D sering dijumpai pada masalah teknik dan fisika. Persamaan differensial Poisson dan Laplace 2D dijumpai pada masalah fluida, potensial, elastisitas, konduksi panas, air tanah dan lain-lain. Seperti persamaan differensial lainnya, kerumitan penyelesaian persamaan differensial Poisson dan Laplace 2D terletak pada bentuk syarat batas yang menyertai persamaan differensial tersebut. Pada tulisan ini diberikan cara penyelesaian persamaan differensial Poisson dan Laplace 2D dengan ketiga tipe syarat batasnya yaitu Dirichlet, Neumann, dan Robbin yang diimplementasikan dalam program komputer menggunakan Matlab.

Implementasi Matlab untuk Menyelesaikan Masalah Syarat Batas Persamaan Differensial Poisson dan Laplace 2d

152

2. Skema Numerik dan Program Komputer a. Partisi Domain Diberikan persamaan differensial Poisson dan Laplace 2D

∂ 2 u ( x, y ) ∂ 2 u ( x, y ) = f ( x, y ) + ∂y 2 ∂x 2

(1)

dengan domain Ω = {( x, y ) | 0 ≤ x ≤ p, 0 ≤ y ≤ q}. Batas-batas domain dapat bertipe Dirichlet, Neumann, atau Robin. Penyelesaian persamaan (1) dengan metode beda hingga, dimulai dengan mempartisi domain seperti pada Gambar 1.

Gambar 1. Domain dan partisinya.

Program komputer (bagian 1): close all; clear all; clc; %Bagian 1 % Menyelesaikan PD Poisson & laplace % dengan domain segi empat dengan 3 tipe syarat batas dapat dipilih % beripe Dirichlet, Neumman, Robbin. % % j


153

% y | Syarat batas 2 % J tinggi ---------------------# % . | | % . | | % . | | % . | | % .Syarat batas 3 | Uxx + Uyy = f | Sayarat batas 4 % . | |% . | | % 3 | | % 2 | | % 1 #---------------------+--- x % 0 lebar % Syarat batas 1 % 1 2 3 ................ I i %-----------------------------------------------------------------disp(' BENTUK UMUM : Uxx + Uyy = f '); disp(' -------------------------------'); disp(' Jenis Persamaan Differensial: '); disp(' 1. Laplace ') disp(' 2. Poisson ') PD=input(' * pilih 1-2 : ? '); if PD==1, f=0; else f=input(' f = ? '); end; %if disp('UKURAN DOMAIN') disp('--------------') lebar= input(' lebar = ? '); tinggi= input(' tinggi = ? '); I=input('I:indeks maks.pilahan pd arah sb-x. i=1,2,.,I. J=input('J:indeks maks.pilahan pd arah sb-y. j=1,2,.,J.

I = ? '); J = ? ');

h=lebar/(I-1); % lebar kisi k=tinggi/(J-1);% tinggi kisi disp(' '); % program bersambung ke bagian berikutnya

b. Skema Iterasi Pada Batas dan Bagian Dalam Domain Persamaan (1) selanjutnya ditulis menjadi

∂ 2u ∂ 2u ( 2 ) i, j + ( 2 ) i, j = f i, j ∂x ∂y

(2)


154

jika digunakan rumus pendekatan

(

u i −1, j − 2u i , j + u i +1, j ∂ 2u ) ≈ 2 i, j ∂x h2

(

u i , j −1 − 2u i , j + u i , j +1 ∂ 2u ) ≈ 2 i, j k2 ∂y

dan

maka persamaan (2) dapat didekati dengan skema

u i −1, j − 2u i , j + u i +1, j h2

+

u i , j −1 − 2u i , j + u i , j +1 k2

= f i, j

atau

u i −1, j + u i +1, j ui, j =

h

2

+

u i. j −1 + u i , j +1

k2 1 · § 1 2¨ 2 + 2 ¸ k ¹ ©h

− f i, j .

(3)

Persamaan (3) merupakan skema untuk mencari u pada titik-titik grid yang terletak pada bagian dalam domain, jadi berindeks i=2,3,...,I-1 j=2,3,...,J-1. Adapun untuk titiktitik grid yang terletak pada batas domain dalam hal ini u i , j dengan salah satu atau kedua indeksnya adalah i=1, i=I, j=1, j=1, j=J, dibutuhkan modifikasi persamaan (3) sesuai syarat batas yang diberikan (Robin atau Neumann). Pada batas bertipe Dirichlet, nilai u

telah

diketahui, sehingga tidak dibutuhkan skema numerik untuk mencari nilai u pada batas tersebut. Pada sudut batas yang dibentuk oleh dua batas bertipe Dirichlet, nilai u diasumsikan sama dengan nilai u rata-ratanya. Jika sudut batas dibentuk oleh batas bertipe Dirichlet dan batas lainnya bertipe Robin atau Neumann maka nilai u pada sudut batas tersebut diasumsikan mengikuti nilai u dari batas Dirichlet. Nilai-nilai u pada batas bertipe Robin atau Neumann belum diketahui, oleh karena itu dibutuhkan skema numerik untuk mencari nilai u pada batas-batas tersebut. Diketahui bentuk umum syarat batas merupakan tipe Robin, yaitu


P u(x,y) + Q

155

∂u ( x, y ) = G ∂n

dengan n adalah vektor arah normal satuan (vektor satuan yang tegak lurus pada kurva batas).

Y

n

n

∂u ∂u = ∂n ∂y

Domain

n

∂u ∂u = ∂n ∂x

∂u ∂u =− ∂n ∂x

n

∂u ∂u =− ∂y ∂n

X

Gambar 2. Vektor normal satuan Jika suatu batas tidak bertipe Dirichlet ( Q ≠ 0 ) maka bentuk syarat batas Robin dapat dinyatakan sebagai

∂u ( x, y ) = α u ( x, y ) + β ∂n

(4)

dengan α = − P / Q , β = G / Q . Jika α = 0 , maka persamaan (4) menjadi syarat batas bertipe Neumann. Skema numerik untuk mencari nilai u pada suatu batas bertipe Robin atau Neumann pada prinsipnya adalah modifikasi persamaan (3) sesuai syarat batas yang diberikan. Hal ini akan dicontohkan pada bagian berikut.


156

∂u = α 2u + β 2 ∂y

Y

q

∂u = α 3u + β 3 ∂x

∂u = α 4u + β 4 ∂x

X

0

∂u = α 1u + β 1 ∂y

p

Gambar 3. Domain dengan syarat batas Robin atau Neumann

Misalkan domain pada Gambar 3.4 dengan batas batas kiri bertipe Robin, yaitu

∂u ( x, y ) = α 3u ( x, y ) + β 3 , x = 0, 0 ≤ y ≤ q ∂x

(5)

dengan α 3 , β 3 bilangan konstan. Persamaan (5) dapat ditulis menjadi

(

∂u ) i , j = α 3 u i , j + β 3 , i = 1, j = 1,2,..., J . ∂x

(6)

Jika digunakan rumus pendekatan

(

u i +1, j − u i −1, j ∂u ) i, j ≈ 2h ∂x

maka persamaan (6) dapat didekati dengan skema

u i +1, j − u i −1, j 2h

= α 3 u i , j + β 3 ; i=1; j=1,2,...,J.

(7)


157

Jika i=1 digunakan pada skema PD Poisson dan Laplace 2D (3) dan persamaan (7) maka akan dijumpai u 0, j dengan j=1,2,...,J. Sedangkan diketahui bahwa indeks terkecil untuk i adalah 1, ini berarti u 0, j adalah titik-titik fiktif . Oleh karena itu u0,j tidak dapat digunakan secara langsung. Persamaan (3) untuk i=1, diperoleh

u 0, j + u 2, j h2

u1, j =

+

u1. j −1 + u1, j +1

k2 1 · § 1 2¨ 2 + 2 ¸ k ¹ ©h

− f1, j .

(8)

Persamaan (7) untuk i=1, diperoleh

u 2, j − u 0, j 2h

= α 3u1, j + β 3

atau

u 0, j = u 2, j − 2hβ 3 − 2hα 3 u1, j .

(9)

Persamaan (9) disubstitusi pada persamaan (8), diperoleh

2u 2, j − 2hβ 3 u1, j =

+

u1. j −1 + u1, j +1

h2 k2 1 · 1 § 1 2¨ 2 + 2 ¸ + 2 2hα 3 k ¹ h ©h

− f1, j .

(10)

Persamaan (10) merupakan skema untuk mencari u1, j , j=2,3,..., J-1. Sedangkan skema untuk mencari u1,1 dan u1, J belum dapat ditentukan, karena melibatkan syarat batas lain yang membentuk sudut-sudut tersebut. Misalkan domain pada Gambar 3.4, batas bagian atas domain juga

bertipe

Robin,yaitu

∂u ( x, y ) = α 2 u ( x, y ) + β 2 , 0 ≤ x ≤ p, y = q ∂y

(11)


158

dengan α 2 , β 2 bilangan konstan. Persamaan (11) dapat ditulis menjadi

(

∂u ) i , j = α 2 u i , j + β 2 , i=1,2,...,I, j=J. ∂y

(12)

Jika digunakan rumus pendekatan

(

u i , j +1 − u i , j −1 ∂u ) i, j ≈ 2k ∂y

maka persamaan (12) dapat didekati dengan

u i , j +1 − u i , j −1 2k

= α 2 u i , j + β 2 ; i=1,2,...,I; j=J

atau

u i ,J +1 = u i ,J −1 + 2kβ 2 + α 2 , i=1,2,...,I .

(13)

Diketahui indeks terbesar untuk j adalah J. Jadi u i , J +1 , i=1,2,3,...,I adalah titik-titik fiktif. Skema PD Poisson (3) untuk j = J adalah

u i −1,J + u i +1,J u i ,J =

h

2

+

u i.J −1 + u i ,J +1

k2 1 · § 1 2¨ 2 + 2 ¸ k ¹ ©h

− f i ,J (14)

Persamaan (13) disubstitusi ke persamaan (14) diperoleh skema untuk u i , J yaitu

u i −1, J − u i +1, J

2u i .J −1 + 2 kβ 2 − f i,J h k2 1 · 1 § 1 2¨ 2 + 2 ¸ − 2 2 k α 2 k ¹ k ©h 2

u i ,J =

+

(15)

dengan i=2,3,...,I-1. Untuk u1,J dan u I,J ( u pada sudut batas domain) masih memerlukan informasi tambahan dari syarat batas yang lain yang membentuk sudut-sudut tersebut. Karena diketahui domain pada Gambar 4, batas kiri dan batas kanan domain bertipe Robin,


159

maka skema u pada sudut kiri atas ( u1,J ) diperoleh dengan mensubstitusi persamaan (9) dan (13) ke persamaan (3), diperoleh

2u 2,J − 2hβ 3

2u i.J −1 + 2kβ 2 − f 1,J 2 h k = . 1 · 1 1 § 1 2¨ 2 + 2 ¸ + 2 2hα 3 − 2 2kα 2 k ¹ h k ©h 2

u1,J

+

(16)

Khusus untuk sudut batas yang dibentuk oleh batas-batas Dirichlet, nilai u pada sudut diasumsikan sama dengan nilai u rata-ratanya. Jika sudut batas dibentuk oleh batas bertipe Dirichlet dan batas yang lainnya bertipe Robin atau Neumann maka nilai u pada sudut tersebut diasumsikan mengikuti nilai u dari batas bertipe Dirichlet. Setelah semua skema untuk u tersedia, maka untuk u yang belum diketahui nilainya diberikan sebarang nilai awal u ( 0 ) , selanjutnya diiterasi menggunakan skema SOR, yaitu

(u i , j ) v = ω.u i , j + (1 − ω ).(u i , j ) v −1

(17)

dengan v nomor iterasi (v = 1,2,3,...), u i , j adalah ui,j dari skema titik grid pada bagian dalam domain (persamaan (3)) dan skema batas Robin atau Neumann. Parameter SOR ω dipilih 0 < ω < 2 . Proses iterasi terus dilakukan sampai memenuhi kriteria konvergensi yang diberikan, yaitu

(u ) − (u )

v −1

v

i, j

dengan

i, j

< ε , ∀i, j

nilai toleransi kekonvergenan biasanya dipilih nilai 0 < ε < 1 .

ketaksamaan (18) dipenuhi maka proses iterasi dihentikan.

Program komputer (bagian 2): % Bagian 2 (sambungan bagian 1) % s : indeks sisi s=1,2,,3,4 disp('SYARAT BATAS'); disp('------------'); for s=1:4

(18) Setelah


fprintf('\n Syarat batas sisi ke');fprintf('%2d ',(s)); fprintf('\n 1. Dirichlet '); fprintf('\n 2. Neumann '); fprintf('\n 3. Campuran (Robyn)\n'); tipe = input(' * pilih 1-3 : '); switch tipe case 1, TipeBC(s)=1; % syarat batas sisi ke s bertipe Dirichlet BC(s)=input(' u = '); case 2, TipeBC(s)=3; % syarat batas bertipe Neummann dianggap alfa(s)=0; % syarat batas bertipe Robyn dgn alfa=0 switch s case 1, du/dy '); fprintf(' case 2, fprintf(' du/dy '); case 3 du/dx '); fprintf(' case 4, du/dx '); fprintf(' otherwise % tidak ada lagi sisi end %switch s BC(s)=input(' = '); %syarat batas sisi ke s bertipe Robyn otherwise, TipeBC(s)=3; switch s case 1, disp(' dU/dy=alfa.u+betta'); case 2, disp(' dU/dy=alfa.u+betta'); case 3 dU/dx=alfa.u+betta'); disp(' case 4, dU/dx=alfa.u+betta'); disp(' otherwise % tidak ada lagi sisi Robyn end %switch s alfa(s)=input(' BC(s)=input(' end % switch tipe end; %for s

alfa ='); betta = ');

160


% inisialisasi U pada batas for i=1:I U(i,1)=0; U(i,J)=0; end; for j=1:J U(1,j)=0; U(I,j)=0; end; % set nilai U pada syarat batas bertipe Dirichlet if TipeBC(1)==1, for i=1:I U(i,1)=BC(1); end; A1=U(1,1); B1=U(I,1); end; % if if TipeBC(2)==1, for i=1:I U(i,J)=BC(2); end; C2=U(1,J); D2=U(I,J); end; % if if TipeBC(3)==1, for j=1:J U(1,j)=BC(3); end; A3=U(1,1); C3=U(1,J); end; % if if TipeBC(4)==1, for j=1:J U(I,j)=BC(4); end; B4=U(I,1); D4=U(I,J); end; % if %----- nilai U pada titik pojok mengikuti sisi Dirichlet----% pojok A if TipeBC(1)==1 & TipeBC(3)==1 , U(1,1)=(A1+A3)/2; end; if TipeBC(1)==1 & TipeBC(3)==3,

161


U(1,1)=A1; end; if TipeBC(1)==3 & TipeBC(3)==1, U(1,1)=A3; end; % pojok B if TipeBC(1)==1 & TipeBC(4)==1, U(I,1)=(B1+B4)/2; end; if TipeBC(1)==1 & TipeBC(4)==3, U(I,1)=B1; end; if TipeBC(1)==3 & TipeBC(4)==1, U(I,1)=B4; end; % pojok C if TipeBC(2)==1 & TipeBC(3)==1, U(1,J)=(C2+C3)/2; end; if TipeBC(2)==1 & TipeBC(3)==3, U(1,J)=C2; end; if TipeBC(2)==3 & TipeBC(3)==1, U(1,J)=C3; end; % pojok D if TipeBC(2)==1 & TipeBC(4)==1, U(I,J)=(D2+D4)/2; end; if TipeBC(2)==1 & TipeBC(4)==3, U(I,J)=D2; end; if TipeBC(2)==3 & TipeBC(4)==1, U(I,J)=D4; end; % set tebakan awal U pada sisi dgn tipe Neumann atau Robyn % sisi ke 1 if TipeBC(1)==3 & TipeBC(2)==1, for i=2:I-1 U(i,1)=(BC(2)-tinggi*BC(1)-tinggi)/(1+tinggi*alfa(1)); end; end; % if

162


% sisi ke 2 if TipeBC(1)==1 & TipeBC(2)==3, for i=2:I-1 U(i,J)=BC(1)+tinggi*BC(2)/(1-tinggi*alfa(2)); end; end; % if % sisi ke 3 if TipeBC(3)==3 & TipeBC(4)==1, for j=2:J-1 U(1,j)=BC(4)-lebar*BC(3)/(1+lebar*alfa(3)); end; end; % if % sisi ke 4 if TipeBC(3)==1 & TipeBC(4)==3, for j=2:J-1 U(I,j)=BC(3)+lebar*BC(4)/(1-lebar*alfa(4)); end; end; % if % Tebakan awal pojok yang dibentuk 2 sisi Neumann atu Robyn % pojok A if TipeBC(1)==3 & TipeBC(3)==3, A1=(U(1,2)-k*BC(1))/(1+k*alfa(1)); A3=U(2,1)-h*BC(3)/(1+h*alfa(3)); U(1,1)=(A1+A3)/2; end; % tebakan awal pojok B if TipeBC(1)==3 & TipeBC(4)==3, B1=U(I,2)-k*BC(1)/(1+k*alfa(1)); B4=U(I-1,1)-h*BC(4)/(1+h*alfa(4)); U(I,1)=(B1+B4)/2; end; % tebakan awal pojok C if TipeBC(2)==3 & TipeBC(3)==3, C2=U(1,J-1)+k*BC(2)/(1-k*alfa(2)); C3=U(2,J)-h*BC(3)/(1+h*alfa(3)); U(1,J)=(C2+C3)/2; end; % tebakan awal pojok D if TipeBC(2)==3 & TipeBC(4)==3, D2=U(I,J-1)+k*BC(2)/(1-k*alfa(2)); D4=U(I-1,J)+h*BC(4)/(1-h*alfa(4)); U(I,j)=(D2+D4)/2;

163


164

end; H=lebar/2; K=tinggi/2; r=sqrt(H*H+K*K); aver = (U(1,1)+U(I,1)+U(1,J)+U(I,J)-r*r*f)/4; %tebakan awal titik interior for i=2:I-1 for j=2:J-1 U(i,j)=0;%aver end end

Ukonv=0; % inisialisasi banyaknya titik konvergen banyakU=I*J; % banyaknya U seluruhnya epsilon=10.^(-6); % nilai toleransi perbedaan nilai U Mj=(cos(pi/(I+1)) + cos(pi/(J+1)))/2; w= 1.5;%2/(1+sqrt(1-Mj.^2)); dh=1/(h*h); dk=1/(k*k); dhk= 2*(dh+dk); iter=0; %-------------------Bagian iterasi-------------------------while Ukonv1&i1&j<J, U(i,j)=(dh*(U(i-1,j)+U(i+1,j))+dk*(U(i,j-1)+U(i,j+1))f)/dhk; %--- skema sisi Neumann atau Robyn (tak termasuk ujungnya)----% sisi ke 1 (sisi bawah) elseif i>1&i1&i1&j<J & TipeBC(3)==3, %sisi ke 3 (sisi kiri)


165

U(1,j)= (dh*(2*U(2,j) - 2*h*BC(3))+dk*(U(1,j1)+U(1,j+1))- f)/(dhk+dh*2*h*alfa(3)); elseif i==I & j>1&j<J & TipeBC(4)==3, %sisi ke 4 (sisi kanan) U(I,j)= (dh*(2*U(I-1,j)+ 2*h*BC(4))+dk*(U(I,j1)+U(I,j+1))- f)/(dhk-dh*2*h*alfa(4)); %-------skema pojok yang dibentuk 2 sisi Robyn----elseif i==1&j==1 & TipeBC(1)==3 & TipeBC(3)==3, % pojok kiri bawah (pojok A) U(1,1)=(dh*(2*U(2,1)-2*h*BC(3))+dk*(2*U(1,2)2*k*BC(1))- f)/... (dhk+dh*2*h*alfa(3)+dk*2*k*alfa(1)); elseif i==I&j==1& TipeBC(1)==3 & TipeBC(4)==3, % pojok kanan bawah (pojok B) U(I,1)=(dh*(2*U(I-1,1)+2*h*BC(4))+dk*(2*U(I,2)2*k*BC(1))- f)/... (dhk-dh*2*h*alfa(4)+dk*2*k*alfa(1)); elseif i==1&j==J & TipeBC(2)==3 & TipeBC(3)==3, % pojok kiri atas (pojok C) U(1,J)=(dh*(2*U(2,J)-2*h*BC(3))+dk*(2*U(1,J1)+2*k*BC(2))- f)/... (dhk+dh*2*h*alfa(3)-dk*2*k*alfa(2)); elseif i==I&j==J & TipeBC(2)==3 & TipeBC(4)==3, % pojok kanan atas (pojok D) U(I,J)=(dh*(2*U(I-1,1)+2*h*BC(4))+dk*(2*U(I,J1)+2*k*BC(2))- f)/... (dhk-dh*2*h*alfa(4)-dk*2*k*alfa(2)); else %--sisi Dirichlet tidak punya skema iterasi %--karena nilai U nya telah diketahui end % if % skema SOR U(i,j)=w*U(i,j)+(1-w)*temp; % kriteria konvergensi beda=abs(U(i,j)-temp); if beda < epsilon, Ukonv=Ukonv+1 ; end;


166

end % for j end% for i fprintf('\nIterasi ke'); fprintf('%5d ',iter); fprintf('banyak U konvergen');fprintf('%5d ',Ukonv); fprintf(' //');fprintf('%5d',banyakU); if iter==100000, break; end %if end % while--------------------End iterasi----------------------% plot nilai U pada domain arsiran= U(1:I,1:J); pcolor(arsiran) colorbar vert shading interp title('Grafik Kontur Solusi Numerik'); xlabel('x'); ylabel('y'); drawnow; fprintf('\nSelesai\n');

3. Simulasi

Diberikan PD Laplace 2D ∂ 2 u ( x , y ) ∂ 2 u ( x, y ) + =0 ∂x 2 ∂y 2 dengan domain Ω = {( x, y ) | 0 ≤ x ≤ 1,0 ≤ y ≤ 1} dan syarat batas:

§ ∂u · § ∂u · § ∂u · u(x,0)= 1, ¨¨ ¸¸ = 2 , ¨ ¸ = 3 u + 1 , ¨ ¸ = −4 u + 1 © ∂x ¹ 0, y © ∂x ¹1, y © ∂y ¹ x ,1 Hasil program komputer diatas : BENTUK UMUM : Uxx + Uyy = f ------------------------------Jenis Persamaan Differensial: 1. Laplace 2. Poisson * pilih 1-2 : ? 1


167

UKURAN DOMAIN -------------lebar = ? 1 tinggi = ? 1 I:indeks maks.pilahan pd arah sb-x. i=1,2,.,I. J:indeks maks.pilahan pd arah sb-y. j=1,2,.,J.

I = ? 33 J = ? 33

SYARAT BATAS -----------Syarat batas sisi ke 1 1. Dirichlet 2. Neumann 3. Campuran (Robyn) * pilih 1-3 : 1 u = 1 Syarat batas sisi ke 2 1. Dirichlet 2. Neumann 3. Campuran (Robyn) * pilih 1-3 : 2 du/dy = 2 Syarat batas sisi ke 3 1. Dirichlet 2. Neumann 3. Campuran (Robyn) * pilih 1-3 : 3 dU/dx=alfa.u+betta alfa =3 betta = 1 Syarat batas sisi ke 4 1. Dirichlet 2. Neumann 3. Campuran (Robyn) * pilih 1-3 : 3 dU/dx=alfa.u+betta alfa =-4 betta = 1 Iterasi Iterasi ... Iterasi Iterasi Iterasi Iterasi Selesai

ke ke

1 1

banyak U konvergen banyak U konvergen

ke ke ke ke

892 893 894 895

banyak banyak banyak banyak

U U U U

konvergen konvergen konvergen konvergen

33 33 1047 1062 1078 1089

// 1089 // 1089 // // // //

1089 1089 1089 1089


168

»

Gambar 4. Solusi numerik Kebenaran solusi numerik diperlihatkan pada Gambar 5 yang menunjukan bahwa solusi numerik sesuai dengan skema numerik dan syarat batas yang diberikan.

Gambar 5. Kebenaran solusi numerik


169

DAFTAR PUSTAKA

[1] Bailey,W. 2003. The SOR algorithm & its Application to Numerical Solution of Eliptic Partial Differential Equation. Ireland: Dublin Institute of Technology. [2] Bassaruddin, T. 1994. Metode Beda Hingga untuk Persamaan Differensial. Jakarta: Elex Media Komputindo. [3] Constantinides, A. 1987. Applied Numerical Methods with Personal Computer. New York:Mc Graw Hill Inc. [4] Nakamura,S. 1991. Applied Numerical Methods with Software. New York: Prentice Hall Inc. [5] Sturler E. 2003. Iterative Methods and Multigrid. http://www.cse.uiuc.edu/cs550/lectures.htm


170

IMPLEMENTASI MATLAB UNTUK MENYELESAIKAN MASALAH SYARAT BATAS PERSAMAAN DIFFERENSIAL POISSON DAN LAPLACE 2D. La Ode Muhammad Umar RRR 1)

Recommend Documents