Tampilkan postingan dengan label Numerical Method. Tampilkan semua postingan
Tampilkan postingan dengan label Numerical Method. Tampilkan semua postingan

Sabtu, 16 Maret 2013

Metode Runge-Kutta Orde 4 Menggunakan MATLAB


Berikut dibahas tentang persamaan diferensial biasa, ordinary differential equations (ODE) yang diklasifikasikan kedalam masalah nilai awal (initial value) dan masalah nilai batas (boundary value), dimana kedua keadaan ini solusinya dispesifikasi pada waktu awal (initial time). Banyak hukum-hukum fisika yang ‘sangat pas ’ diformulasikan dalam bentuk persamaan diferensial. Lebih lanjut, tidak mengherankan bahwa solusi komputasi numerik dari persamaan-persamaan diferensial menjadi bagian yang umum dalam pemodelan sistem-sistem fisika. Beberapa hukum mendasar diantaranya sebagai berikut:


METODE EULER

Sistem persamaan diferensial yang kompleks dapat dengan mudah dipecahkan menggunakan metode numerik, salah satu metode yang paling sederhana adalah metode Euler.

Pendekatan metode Euler menggunakan pendekatan ekspansi Taylor berikut:


Dengan orde yang besar dari persamaan diatas nilainya kecil, maka Euler hanya menggunakan sampai orde 1 saja, sehingga diperoleh


dengan h adalah step size.

Metode Euler ini memiliki kelemahan, untuk t yang semakin besar, nilai eror dari solusi akan semakin besar, untuk itu lahirlah metode Runge-Kutta yang menawarkan penyelesaian persamaan diferensial dengan pertumbuhan truncation error yang jauh lebih kecil.

METODE RUNGE KUTTA ORDE 4

Persamaan Runge-Kutta orde 4 dapat dituliskan sebagai berikut:


CONTOH KASUS SISTEM DIFERENSIAL BIASA

Contoh aplikasi sistem persamaan diferensial orde satu terkopel adalah persamaan Lorenz tahun 1963, yang menggambarkan fenomena konveksi udara yang dibangun atas tiga persamaan diferensial terkopel berikut

dimana a adalah bilangan Prandtl, r adalah bilangan Rayleighr, dan b adalah geometris dari sistem fisik yang ditinjau. Sedangkan x adalah laju aliran konveksi. y adalah perbedaan temperature horizontal aliran konveksi dan z perbedaan temperature horizontal aliran konveksi terhadap titik equilibrium.

Untuk dapat memahami fenomena dari persamaan Lorenz diatas, maka persamaan tersebut harus dicari solusinya. Berikut adalah solusi untuk sistem persamaan Lorenz diatas menggunakan MATLAB.

Metode Numerik Sistem Lorenz
Program berikut dibuat dalam tiga file terpisah, yaitu file fungsi, file metode numerik, dan file untuk memplot grafik karena pada MATLAB memiliki fasilitas dapat memanggil file. Sebenarnya program juga dapat dibuat hanya dalam satu file saja, tetapi ketika kita akan membuat fungsi baru program menjadi tidak efisien karena semuanya harus diubah, tetapi dengan sistem panggil file, ketika kita akan merubah fungsi atau output tidak perlu lagi untuk merubah metode numeriknya.

Fungsi M-File Sistem Lorenz

function f = lorenz(t,y)

f = zeros(1,3);
f(1) = 10 * (y(2) - y(1));
f(2) = -y(1) * y(3) + 28 * y(1) - y(2);
f(3) = y(1) * y(2) - 8 * y(3) / 3;


Fungsi Metode Runge Kutta  Orde 4
function [tSol,ySol] = RK4(dEqs,t,y,tStop,h)
% Cara Menggunakan Metode Runge Kutta Orde Empat
% Cara Menggunakan: [tSol,ySol] = RK4(@namafungsidiferensial,t,y,tStop,h)
% INPUT:
% namafungsidiferensial= merupakan sistem pers.diferensial orde satu
% F(t,y) = [dy1/dt dy2/dt dy2/dt ...].
% t,y = kondisi awal; y harus dibuat dalam vektor kolom.
% tStop = batas akhir waktu
% h = step size
% OUTPUT:
% t = variabel bebas (waktu)
% y = variabel yang bergantung waktu
if size(y,1) > 1 ;
    y = y';
end
tSol = zeros(2,1); ySol = zeros(2,length(y));
tSol(1) = t; ySol(1,:) = y;
i = 1;
while t < tStop
i = i + 1;
h = min(h,tStop - t);
K1 = h*feval(dEqs,t,y);
K2 = h*feval(dEqs,t + h/2,y + K1/2);
K3 = h*feval(dEqs,t + h/2,y + K2/2);
K4 = h*feval(dEqs,t + h,y + K3);
y = y + (K1 + 2*K2 + 2*K3 + K4)/6;
t = t + h;
tSol(i) = t; ySol(i,:) = y; % Untuk mengeluarkan solusi
end

Script Eksekusi Program M-File Sistem Lorenz
[t,y] = RK4(@lorenz,0,[0.1 0.1 0.1],50,0.01);
figure(1)
plot3(y(:,1),y(:,2),y(:,3))
xlabel('x')
ylabel('y')
zlabel('z')
title('DIAGRAM FASE SISTEM KONVEKSI LORENZ')
figure (2)
plot(t,y)
xlabel('Waktu')
ylabel('Dinamika')
title('TIME SERIES SISTEM KONVEKSI LORENZ')
legend('x','y','z')

Output Program

Gambar 1. Diagram fase sistem Lorenz

Gambar 2. Grafik time series sistem Lorenz

Dari gambar 1 dan 2 terlihat bahwa sistem Lorenz memiliki solusi yang tidak periodik, solusi ini dikenal dengan fenomena CHAOS yang memiliki ciri sulit untuk diprediksi untuk interval waktu yang lama. dan ketika kita mengubah sedikit saja salah satu parameter atau kondisi awal sistem, maka akan melahirkan solusi yang sangat berbeda untuk selang waktu yang lama.

Copyright 2013 @ Profesor Bolabot














Share:
Read More