Listing Program C untuk Heat Equation 1-D dengan Syarat Batas Neumann


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/2 satuan 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:
          Teknik Iterasi
          Install Gnuplot di Ubuntu
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.

Solusi Numerik Finite Difference Skema Eksplisit
Baca juga: Diskritisasi Heat Equation 2-D (Skema Eksplisit)
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 step k + 1 secara eksplisit karena nilai-nilai u di sepanjang batang pada time step k  (terutama pada k = 0) sudah diketahui dari syarat awal.
• Diskritisasi syarat awal
• Diskritisasi syarat batas (syarat batas Neumann)
• Syarat stabilitas
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)
Bandingkan dengan grafik solusi analitik.

Solusi Numerik Finite Difference Skema Implisit
Baca juga: Diskritisasi Persamaan Laplace 2-D (Skema Implisit)
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 step k + 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.

Solusi Numerik Finite Difference Skema Crank-Nicolson
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.
          BACA JUGA:



2 comments :

  1. Terima kasih banyak..
    Sangat membantu..

    ReplyDelete
  2. Kereen...
    Izin seraapp..

    ReplyDelete

Silahkan Tulis Komentar Kamu :)