Sebenarnya artikel kali ini adalah pesanan dari beberapa
orang. Masalah yang dihadapi adalah bagaimana cara memecahkan heat equation pada batang 1-D dengan
syarat batas Neumann di kedua ujung-ujungnya. Heat equation dipecahkan secara analitik dan numerik, kemudian
dibandingkan antara solusi analitik dengan numeriknya untuk interval (Δx) yang sama. Solusi numerik dihitung
menggunakan metode finite difference
dengan skema eksplisit, implisit, dan Crank-Nicolson, dengan bantuan komputer
tentunya.
Adapun masalahnya adalah sebagai berikut
Heat equation
Syarat batas Neumann untuk batang sepanjang L
Syarat awal (u di sepanjang batang pada t
= 0)
Sebenarnya, memecahkan masalah di atas tidaklah sulit,
dengan pertimbangan bahwa kita hanya memecahkan heat equation 1-dimensi dengan syarat batas Neumann dan syarat awal
seperti di atas. Kita anggap u(x,t)
adalah temperatur.
Jika kita lihat syarat awalnya, dengan mudah kita dapat
melihat bahwa pada saat awal (t =0), grafik temperatur di sepanjang
batang adalah seperti berikut ini
Gambar 1. Grafik Temperatur di Sepanjang Batang Saat t = 0
Secara intuisi, bisa dilihat pula bahwa temperatur u di sepanjang batang semakin lama akan
menuju L/2satuan temperatur, karena syarat batas Neumann pada kedua ujung batang.
Akan kita lihat pada pembahasan selanjutnya bahwa intuisi ini benar untuk kasus
heat equation di sini.
Sebelum lebih jauh ke pembahasan solusi dan listing
program C untuk menghitung solusi heat
equation di atas, sebaiknya baca dulu artikel berikut ini untuk pendahuluan:
Artikel di atas penting karena kita akan membahas tentang
iterasi persamaan linier dalam bentuk matriks untuk menghitung solusi heat equation di atas (detil-detil
perhitungannya ada di artikel-artikel tersebut); kita juga akan menggunakan
listing program C yang digabungkan dengan pipe
Gnuplot agar hasil hitung kita dapat langsung secara otomatis diplot ke dalam
bentuk grafik setelah listing program C tersebut dikompilasi dan dijalankan.
Penting: Sebaiknya gunakan
Ubuntu (atau OS berbasis UNIX lainnya) yang sudah diinstalasi Gnuplot untuk
mengkompilasi dan menjalankan listing program C di artikel ini karena listing
program C tersebut kami tulis dan jalankan di Ubuntu pula.
Jika pembaca sudah familiar dengan teknik iterasi dan pipe Gnuplot pada bahasa pemrograman C, selanjutnya
kita akan membahas solusi heat equation
1-D untuk kasus kita di atas. Di sini kita tidak akan menghitung error solusi numerik terhadap solusi
analitik dalam persentase, yang akan kita lakukan di sini adalah membandingkan
grafik solusi analitik dengan grafik solusi numerik secara visual agar lebih
menarik. Error dalam persentase dapat
dihitung dengan mudah dengan program C di sini, sehingga hal tersebut kami
serahkan kepada pembaca.
Solusi Analitik
Solusi analitik (solusi khusus) heat equation 1-D di atas dapat dibuktikan dengan mudah, yang mana
solusi analitiknya adalah sebagai berikut (tidak perlu diturunkan di sini)
Yang akan kita lakukan selanjutnya adalah memplot solusi
analitik di atas ke dalam bentuk grafik. Untuk memplot solusi analitik di atas,
kita perlu menghitung solusi analitiknya terlebih dahulu. Untuk menghitung dan memplot
solusi analitik di atas, dapat digunakan listing program C berikut ini (compile di Ubuntu yang sudah diintegrasikan
dengan gnuplot; compile dengan gcc
standar):
/* FISIKA VERITAS
fisikaveritas.blogspot.com
Heat equation 1-D, Syarat Batas Neumann
Plot Solusi Analitik */
#include "stdio.h"
#include "stdlib.h"
#include "math.h"
main()
{
int n, i, k, N, tmax;
float u[50][50], deret[50][50], deretan[50][50], L, dt, alpha, pi, dx;
FILE *simpan;
simpan = fopen("data-an.dat", "w");
pi = 3.1416;
/* Bagian yang dapat diubah-ubah */
L = 10; // Panjang Batang
N = 10;
alpha = 1;
dt = 0.5;
dx = 1;
tmax = 7;
printf("\n\n\t PERHITUNGAN ANALITIK HEAT EQUATION \n\n\n");
/* PERHITUNGAN ANALITIK */
for(k=0; k<=tmax; k++)
{
for(i=0; i<=N; i++)
{
for(n=1; n<=500; n++) // Penjumlahan Sigma, Sampai n = 500;
{
deret[i][k] = (cos((((2*n)-1)*pi*i*dx)/(L))*exp(-(alpha*alpha*((2*n)-1)*((2*n)-1)*pi*pi*k*dt)/(L*L)))/(((2*n)-1)*((2*n)-1));
deretan[i][k] += deret[i][k];
}
u[i][k] = L/2 - deretan[i][k]*((4*L)/(pi*pi));
}
}
/* TAMPILAN PROGRAM DALAM BENTUK DATA */
for(i=0; i<=N; i++)
{
printf(" x = %d; ", i);
for(k=0; k<=tmax; k++)
{
printf(" %.2f ", u[i][k]);
fprintf(simpan," %.2f", u[i][k]);
}
printf("\n");
fprintf(simpan, "\n");
}
printf("\n\n\n");
/* Menampilkan Plot Data dengan Gnuplot */
FILE *pipe = popen("gnuplot -persist","w");
fprintf(pipe, "set title 'Heat Equation Analitik: Grafik Temperatur Sepanjang Batang'\n");
fprintf(pipe, "set xlabel 'x' \n");
fprintf(pipe, "set xrange [0:14.5]\n");
fprintf(pipe, "set ylabel 'Temperatur' \n");
fprintf(pipe, "plot 'data-an.dat' using 0:1 with lines title 't = 0.0 detik', '' using 0:2 with lines title 't = 0.5 detik', '' using 0:3 with lines title 't = 1.0 detik', '' using 0:3 with lines title 't = 1.5 detik', '' using 0:4 with lines title 't = 2.0 detik', '' using 0:5 with lines title 't = 2.5 detik', '' using 0:6 with lines title 't = 3.0 detik' \n");
// gunakan bagian ini jika ingin tampilan plot dalam bentuk kurva dengan titik.
/*
fprintf(pipe, "plot 'data-an.dat' using 0:1 with linespoints title 't = 0.0 detik', '' using 0:2 with linespoints title 't = 0.5 detik', '' using 0:3 with linespoints title 't = 1.0 detik', '' using 0:3 with linespoints title 't = 1.5 detik', '' using 0:4 with linespoints title 't = 2.0 detik', '' using 0:5 with linespoints title 't = 2.5 detik', '' using 0:6 with linespoints title 't = 3.0 detik' \n");
*/
close(pipe);
return 0;
}
Pembaca dapat mengubah-ubah parameter sesuai yang
diinginkan langsung di dalam listing program C tersebut. Untuk perhitungan di
sini, diambil selang waktu 0.5 detik, kita ambil pula deret penjumlahan sigma
di atas sebanyak 500 deret n.
Biasanya, semakin besar n yang
diambil maka hasil hitung menjadi semakin teliti pula. Adapun grafik solusi
analitiknya setelah listing programnya dikompilasi dan dijalankan adalah
sebagai berikut
Gambar 2. Grafik Solusi Analitik (dengan Gnuplot)
Terlihat bahwa pada t
=0, grafik di atas menunjukkan temperatur
di sepanjang pelat adalah seperti pada gambar 1; tentu saja, karena syarat awal
kasus heat equation di atas adalah u(x,0)
= x.
Kita akan menghitung solusi numerik persamaan di atas
dengan menggunakan metode finite
difference dengan skema eksplisit. Jika kita menggunakan skema eksplisit
untuk menghitung solusi persamaan diferensial, maka kita akan berhadapan dengan
syarat stabilitas, maksudnya, kita tidak bisa sembarang memilih selang waktu (time step size) Δt, namun harus memenuhi syarat stabilitas tersebut agar error perhitungan solusi semakin kecil seiring kenaikan time step.
• Kita lakukan diskritisasi heat equation
Kita tulis
Jadi didapatkan matriks heat equation sebagai berikut
Bentuk di atas dikatakan eksplisit karena
kita dapat secara langsung menghitung u
pada time stepk + 1 secara eksplisit karena nilai-nilai u di sepanjang batang pada time step k (terutama pada k = 0) sudah diketahui dari syarat awal.
Sekarang kita kumpulkan semua persamaan
diskrit di atas seperti berikut
Stencil untuk bentuk diskrit heat equation skema eksplisit adalah
sebagai berikut
Gambar 3. Stencil Heat Equation Skema Eksplisit
Adapun listing program C (dengan Gnuplot)
untuk menghitung solusi heat equation
menggunakan skema eksplisit di atas adalah sebagai berikut
/* FISIKA VERITAS
fisikaveritas.blogspot.com
Heat equation 1-D, Syarat Batas Neumann
Finite Difference Skema Eksplisit */
#include "stdio.h"
#include "stdlib.h"
main()
{
int i, k, L, N;
float u[50][50], alpha, g, dt, dx, tmax;
FILE *simpan;
simpan = fopen("data-e.dat", "w");
L = 10; //Panjang batang
N = 10; //Jumlah grid (pilih N kurang dari atau sama dengan L)
alpha = 1;
tmax = 7; // tmax adalah time step maksimum
/* BENTUK DISKRIT SUMBU-X */
dx = L/N;
/* SYARAT STABILITAS */
dt = (dx*dx)/(2*alpha);
/* GAMMA g */
g = ((alpha*dt)/(dx*dx));
/* SYARAT AWAL */
for(i=0; i<=N; i++)
{
u[i][0] = i; // u(x,0) = x
}
/* PERHITUNGAN */
for(k=0; k<=tmax; k++)
{
/* SYARAT BATAS NEUMANN */
u[-1][k] = u[0][k]; // backward difference
u[N+1][k] = u[N][k]; // forward difference
for(i=0; i<=N; i++)
{
/* PERHITUNGAN HEAT EQUATION */
u[i][k+1] = g*(u[i+1][k] - 2*u[i][k] + u[i-1][k]) + u[i][k];
}
}
/* TAMPILAN PROGRAM DALAM BENTUK DATA */
printf("\n\n\n\t HASIL PERHITUNGAN HEAT EQUATION SKEMA EKSPLISIT\n\n");
printf(" t =");
for(k=0; k<=tmax; k++)
{
printf(" %.2f ", (dt*k));
}
printf("\t detik \n -------------------------------------------------------\n");
for(i=0; i<=N; i++)
{
printf(" x = %d; ", i);
for(k=0; k<=tmax; k++)
{
printf(" %.2f ", u[i][k]);
fprintf(simpan," %.2f", u[i][k]);
}
printf("\n");
fprintf(simpan, "\n");
}
printf("\n\n\n");
/* Menampilkan Plot Data dengan Gnuplot */
FILE *pipe = popen("gnuplot -persist","w");
fprintf(pipe, "set title 'Heat Equation Skema Eksplisit: Grafik Temperatur Sepanjang Batang'\n");
fprintf(pipe, "set xlabel 'x' \n");
fprintf(pipe, "set xrange [0:14.5]\n");
fprintf(pipe, "set ylabel 'Temperatur' \n");
fprintf(pipe, "plot 'data-e.dat' using 0:1 with linespoints title 't = 0.0 detik', '' using 0:2 with linespoints title 't = 0.5 detik', '' using 0:3 with linespoints title 't = 1.0 detik', '' using 0:3 with linespoints title 't = 1.5 detik', '' using 0:4 with linespoints title 't = 2.0 detik', '' using 0:5 with linespoints title 't = 2.5 detik', '' using 0:6 with linespoints title 't = 3.0 detik' \n");
// gunakan bagian ini jika ingin tampilan plot dalam bentuk titik-titik.
/*
fprintf(pipe, "plot 'data-e.dat' using 0:1 with points title 't = 0.0 detik', '' using 0:2 with points title 't = 0.5 detik', '' using 0:3 with points title 't = 1.0 detik', '' using 0:3 with points title 't = 1.5 detik', '' using 0:4 with points title 't = 2.0 detik', '' using 0:5 with points title 't = 2.5 detik', '' using 0:6 with points title 't = 3.0 detik' \n");
*/
close(pipe);
return 0;
}
Hasil perhitungan untuk Δx = 1, jumlah grid N =
10, dan Δt = 0.5 detik dengan
menggunakan listing program di atas (setelah dikompilasi dan dijalankan)
ditunjukkan oleh grafik hasil hitung berikut
Gambar 4. Grafik Solusi Numerik Skema Eksplisit (dengan Gnuplot)
Keuntungan menggunakan skema implisit untuk menghitung
solusi numerik suatu persamaan diferensial parsial adalah tidak diperlukannya
syarat stabilitas, jadi kita bisa mengambil berapapun besarnya Δt untuk perhitungan solusi. Namun
kekurangannya dibandingkan dengan skema eksplisit adalah, algoritma pemrograman
menjadi sedikit lebih rumit karena melibatkan iterasi matriks; dan juga perhitungan
solusi lebih lama daripada perhitungan solusi skema eksplisit (untuk time step size yang sama) karena untuk
setiap time step k, kita perlu
mengiterasi matriks persamaan diskrit heat
equation. Error perhitungan pada
skema implisit bergantung pada jumlah iterasi yang kita gunakan; biasanya,
semakin banyak iterasi yang dilakukan maka error
menjadi semakin kecil.
• Kita lakukan diskritisasi heat equation (perbedaan dengan skema eksplisit hanya pada
pengambilan time step untuk turunan
kedua terhadap x)
Kita tulis
Jadi didapatakan matriks heat equation sebagai berikut
Bentuk di atas dikatakan implisit karena kita
tidak dapat menghitung secara langsung u
pada time stepk + 1, karena nilai-nilai u
di sepanjang batang pada time step k +
1 belum diketahui, satu-satunya nilai
yang sudah diketahui hanyalah u pada time step k. Jika kita mengambil jumlah
grid untuk perhitungan sebanyak N,
maka dari matriks heat equation di
atas, ada sebanyak N buah persamaan
dan juga N buah variabel yang akan
dipecahkan; artinya kita bisa menggunakan iterasi dengan tebakan awal
tertentu untuk menghitung solusinya. Di sini kita menggunakan iterasi Gauss-Seidell.
Syarat awal dan syarat batas Neumann sama seperti pada
skema eksplisit.
Stencil untuk bentuk diskrit heat equation skema implisit adalah
sebagai berikut
Gambar 5. Stencil Heat Equation Skema implisit
Adapun listing program C (dengan Gnuplot)
untuk menghitung solusi heat equation
menggunakan skema implisit di atas adalah sebagai berikut
/* FISIKA VERITAS
fisikaveritas.blogspot.com
Heat equation 1-D, Syarat Batas Neumann
Finite Difference Skema Implisit */
#include "stdio.h"
#include "stdlib.h"
main()
{
int i, k, L, N, iterasi;
float u[50][50], alpha, g, dt, dx, tmax;
FILE *simpan;
simpan = fopen("data-im.dat", "w");
L = 10; //Panjang batang
N = 10; //Jumlah grid (pilih N kurang dari atau sama dengan L)
/* BENTUK DISKRIT SUMBU-X */
dx = L/N;
/* ALPHA DAN dt */
alpha = 1;
dt = 0.5;
/* GAMMA g */
g = ((alpha*dt)/(dx*dx));
/* T max */
tmax = 7; // tmax adalah time step maksimum
/* SYARAT AWAL */
for(i=0; i<=N; i++)
{
u[i][0] = i; // u(x,0) = x
}
/* SYARAT BATAS NEUMANN UNTUK T = 0 */
u[-1][0] = u[0][0]; // backward difference
u[N+1][0] = u[N][0]; // forward difference
/* PERHITUNGAN */
for(k=0; k<=tmax; k++)
{
/* TEBAKAN AWAL */
u[i+1][k+1] = 5;
u[i-1][k+1] = 5;
for(iterasi=1; iterasi<=100; iterasi++) // Iterasi GAUSS-SEIDEL
{
/* SYARAT BATAS NEUMANN HARUS TERJAGA DI SETIAP ITERASI */
u[-1][k+1] = u[0][k+1]; // backward difference
u[N+1][k+1] = u[N][k+1]; // forward difference
for(i=0; i<=N; i++)
{
/* MATRIKS PERHITUNGAN HEAT EQUATION */
u[i][k+1] = (g/(1+(2*g)))*(u[i+1][k+1] + u[i-1][k+1]) + (1/(1+(2*g)))*u[i][k];
}
}
}
/* TAMPILAN PROGRAM DALAM BENTUK DATA */
printf("\n\n\n\t HASIL PERHITUNGAN HEAT EQUATION SKEMA IMPLISIT\n\n");
printf(" t =");
for(k=0; k<=tmax; k++)
{
printf(" %.2f ", (dt*k));
}
printf("\t detik \n -------------------------------------------------------\n");
for(i=0; i<=N; i++)
{
printf(" x = %d; ", i);
for(k=0; k<=tmax; k++)
{
printf(" %.2f ", u[i][k]);
fprintf(simpan," %.2f", u[i][k]);
}
printf("\n");
fprintf(simpan, "\n");
}
printf("\n\n\n");
/* Menampilkan Plot Data dengan Gnuplot */
FILE *pipe = popen("gnuplot -persist","w");
fprintf(pipe, "set title 'Heat Equation Skema Implisit: Grafik Temperatur Sepanjang Batang'\n");
fprintf(pipe, "set xlabel 'x' \n");
fprintf(pipe, "set xrange [0:14.5]\n");
fprintf(pipe, "set ylabel 'Temperatur' \n");
fprintf(pipe, "plot 'data-im.dat' using 0:1 with linespoints title 't = 0.0 detik', '' using 0:2 with linespoints title 't = 0.5 detik', '' using 0:3 with linespoints title 't = 1.0 detik', '' using 0:3 with linespoints title 't = 1.5 detik', '' using 0:4 with linespoints title 't = 2.0 detik', '' using 0:5 with linespoints title 't = 2.5 detik', '' using 0:6 with linespoints title 't = 3.0 detik' \n");
// gunakan bagian ini jika ingin tampilan plot dalam bentuk titik-titik.
/*
fprintf(pipe, "plot 'data-e.dat' using 0:1 with points title 't = 0.0 detik', '' using 0:2 with points title 't = 0.5 detik', '' using 0:3 with points title 't = 1.0 detik', '' using 0:3 with points title 't = 1.5 detik', '' using 0:4 with points title 't = 2.0 detik', '' using 0:5 with points title 't = 2.5 detik', '' using 0:6 with points title 't = 3.0 detik' \n");
*/
close(pipe);
return 0;
}
Hasil perhitungan untuk Δx = 1, jumlah grid N =
10, Δt = 0.5 detik, dan iterasi sebanyak 100 kali (untuk setiap time step) dengan
menggunakan listing program di atas (setelah dikompilasi dan dijalankan)
ditunjukkan oleh grafik hasil hitung berikut
Gambar 6. Grafik Solusi Numerik Skema Implisit (dengan Gnuplot)
Bandingkan dengan grafik solusi analitik dan solusi numerik
skema eksplisit.
Skema Crank-Nicolson sebetulnya adalah perata-rataan
skema Implisit dengan skema Eksplisit, bentuknya untuk heat equation adalah sebagai berikut
Bisa juga kita tulis lebih umum seperti berikut ini
Jika
θ = 0 ; maka bentuk di
atas menjadi skema Implisit
θ = 1 ; maka bentuk di
atas menjadi skema Eksplisit
θ = ½ ; maka bentuk di
atas menjadi skema Crank-Nicolson
Kita tulis lagi bentuk umumnya sebagai berikut
dengan
Dari persamaan di atas, dapat kita tulis bentuk yang akan
kita gunakan untuk menulis listing program C sebagai berikut
Syarat awal dan syarat batas Neumann sama
seperti sebelumnya. Perhitungan menggunakan iterasi Gauss-Seidell.
Stencil untuk bentuk diskrit heat equation skema Crank-Nicolson
adalah sebagai berikut
Gambar 7. Stencil Heat Equation Skema Crank-Nicolson
Adapun listing program C (dengan Gnuplot)
untuk menghitung solusi heat equation
menggunakan skema Crank-Nicolson di atas adalah sebagai berikut
/* FISIKA VERITAS
fisikaveritas.blogspot.com
Heat equation 1-D, Syarat Batas Neumann
Finite Difference Skema Crank-Nicolson */
#include "stdio.h"
#include "stdlib.h"
main()
{
int i, k, L, N, iterasi;
float u[50][50], alpha, g, dt, dx, tmax, theta;
FILE *simpan;
simpan = fopen("data-CN.dat", "w");
L = 10; //Panjang batang
N = 10; //Jumlah grid (pilih N kurang dari atau sama dengan L)
/* BENTUK DISKRIT SUMBU-X */
dx = L/N;
/* THETA
rentangnya adalah dari 0 sampai 1.
jika theta = 0, maka bentuk diskrit kembali menjadi skema implisit,
jika theta = 1, maka bentuk diskrit kembali menjadi skema eksplisit,
pada Skema Crank-Nicolson ini kita ambil theta = 0.5 */
theta = 0.5;
/* ALPHA, dt */
alpha = 1;
// dt = (dx*dx)/(2*alpha); /* dipakai jika theta diambil sama dengan 1 (Skema eksplisit) */
dt = 0.5;
/* GAMMA g */
g = ((alpha*dt)/(dx*dx));
/* T max */
tmax = 7; // tmax adalah time step maksimum
/* SYARAT AWAL */
for(i=0; i<=N; i++)
{
u[i][0] = i; // u(x,0) = x
}
/* SYARAT BATAS NEUMANN UNTUK T = 0 */
u[-1][0] = u[0][0]; // backward difference
u[N+1][0] = u[N][0]; // forward difference
/* PERHITUNGAN */
for(k=0; k<=tmax; k++)
{
/* TEBAKAN AWAL */
u[i+1][k+1] = 5;
u[i-1][k+1] = 5;
for(iterasi=1; iterasi<=100; iterasi++) // Iterasi GAUSS-SEIDEL
{
/* SYARAT BATAS NEUMANN HARUS TERJAGA DI SETIAP ITERASI */
u[-1][k+1] = u[0][k+1]; // backward difference
u[N+1][k+1] = u[N][k+1]; // forward difference
for(i=0; i<=N; i++)
{
/* MATRIKS PERHITUNGAN HEAT EQUATION */
u[i][k+1] = (u[i][k]*(1-(2*g*theta)) + (u[i+1][k+1] + u[i-1][k+1])*(g*(1-theta)) + (u[i+1][k]+u[i-1][k])*(g*theta)) / (1+(2*g*(1-theta)));
}
}
}
/* TAMPILAN PROGRAM DALAM BENTUK DATA */
printf("\n\n\n\t HASIL PERHITUNGAN HEAT EQUATION SKEMA CRANK-NICOLSON\n\n");
printf(" t =");
for(k=0; k<=tmax; k++)
{
printf(" %.2f ", (dt*k));
}
printf("\t detik \n -------------------------------------------------------\n");
for(i=0; i<=N; i++)
{
printf(" x = %d; ", i);
for(k=0; k<=tmax; k++)
{
printf(" %.2f ", u[i][k]);
fprintf(simpan," %.2f", u[i][k]);
}
printf("\n");
fprintf(simpan, "\n");
}
printf("\n\n\n");
/* Menampilkan Plot Data dengan Gnuplot */
FILE *pipe = popen("gnuplot -persist","w");
fprintf(pipe, "set title 'Heat Equation Skema Crank-Nicolson: Grafik Temperatur Sepanjang Batang'\n");
fprintf(pipe, "set xlabel 'x' \n");
fprintf(pipe, "set xrange [0:14.5]\n");
fprintf(pipe, "set ylabel 'Temperatur' \n");
fprintf(pipe, "plot 'data-CN.dat' using 0:1 with linespoints title 't = 0.0 detik', '' using 0:2 with linespoints title 't = 0.5 detik', '' using 0:3 with linespoints title 't = 1.0 detik', '' using 0:3 with linespoints title 't = 1.5 detik', '' using 0:4 with linespoints title 't = 2.0 detik', '' using 0:5 with linespoints title 't = 2.5 detik', '' using 0:6 with linespoints title 't = 3.0 detik' \n");
// gunakan bagian ini jika ingin tampilan plot dalam bentuk titik-titik.
/*
fprintf(pipe, "plot 'data-e.dat' using 0:1 with points title 't = 0.0 detik', '' using 0:2 with points title 't = 0.5 detik', '' using 0:3 with points title 't = 1.0 detik', '' using 0:3 with points title 't = 1.5 detik', '' using 0:4 with points title 't = 2.0 detik', '' using 0:5 with points title 't = 2.5 detik', '' using 0:6 with points title 't = 3.0 detik' \n");
*/
close(pipe);
return 0;
}
Hasil perhitungan untuk Δx = 1, jumlah grid N =
10, Δt = 0.5 detik, dan iterasi sebanyak 100 kali (untuk setiap time step) dengan
menggunakan listing program di atas (setelah dikompilasi dan dijalankan)
ditunjukkan oleh grafik hasil hitung berikut
Gambar 8. Grafik Solusi Numerik Skema Crank-Nicolson
Bandingkan dengan grafik solusi analitik, solusi numerik
skema eksplisit, dan skema implisit.
Sekarang kita dapat membandingkan seluruh grafik hasil
perhitungan untuk kasus heat equation
1-D di atas secara visual untuk setiap iΔx yang sama; seluruh grafik hasil hitung
ditunjukkan oleh gambar 9 berikut ini. (Klik gambar untuk memperbesar tampilan gambar)
Gambar 9. Perbandingan Semua Metode Solusi
Untuk justifikasi hasil hitung, kami serahkan kepada
pembaca. Semoga bermanfaat.
Download:
Jika listing program yang tercantum di
artikel ini tidak dapat di-copy atau
bermasalah (karena formatnya berubah sehingga tidak dapat di-compile), silahkan download listing
programnya di sini: heat-1D.rar(5,20
kb)
TRIVIA:
Secara intuitif, bisa kita amati untuk kasus
di atas, temperatur di sepanjang batang akan seragam seiring waktu, yaitu
temperaturnya mendekati L/2. Pada
contoh perhitungan kita di atas, panjang batang adalah 10 satuan panjang,
artinya, temperatur batang seiring waktu akan seragam mendekati 5 satuan
temperatur.
Kita lihat juga dari solusi analitik di atas,
jika t ⟶ ∞, maka suku pada deret sigma akan mendekati 0, yang tersisa adalah suku L/2. Jadi untuk kasus kita
Jika kita plot grafik analitiknya seiring
waktu dengan pengambilan Δt yang
bersesuaian dan untuk jumlah grid N =
10, maka hasilnya adalah sebagai berikut.
Gambar 10. Grafik Solusi Analitik Sampai 30 detik Pertama Menunjukkan Bahwa Batang Hampir Steady State
Dalam bentuk animasi:
Gambar 11. Grafik Solusi Analitik Sampai 25 detik Pertama Menunjukkan Bahwa Batang Hampir Steady State
Jika solusi numerik diplot dengan paramater
yang sama dengan plot solusi analitik di atas (kami serahkan kepada pembaca
untuk mencoba-coba listing programnya), maka hasilnya hampir sama dengan solusi
analitik. Solusi analitik ini digunakan untuk memvalidasi solusi numerik kita,
jadi dapat kita simpulkan bahwa hasil hitung solusi numerik sudah mendekati
kebenaran.
Terima kasih banyak..
ReplyDeleteSangat membantu..
Kereen...
ReplyDeleteIzin seraapp..