Bagian 2: Implementasi Simulasi

Bagian 1 meninggalkan persamaan di atas kertas; sekarang kita jalankan. Dua berkas menjadi tulang punggung tahap ini: modelPertumbuhan.m yang menghitung turunan, dan jalankanSimulasi.m yang memanggil solver ode45 dengan pengaturan yang membuat hasilnya stabil dan bisa diuji. Setiap keputusan implementasi — bentuk argumen, vektor waktu, toleransi — diambil agar hasil simulasi layak dipercaya dan layak dibandingkan dengan solusi analitik.

Yang kita bangun pada bagian ini:

• modelPertumbuhan.m — fungsi ODE tiga argumen,
• jalankanSimulasi.m — panggilan ode45 deterministik dengan odeset,
• validasi terhadap solusi analitik lewat assert,
• struktur hasil berupa struct dengan t, P, dan waktuSetengahK.

Fungsi Model: modelPertumbuhan.m #

Solver ODE di MATLAB tidak menerima persamaan sebagai simbol; solver menerima fungsi yang diberi waktu t dan keadaan P, lalu mengembalikan turunan dP/dt. Berikut seluruh isi src/modelPertumbuhan.m:

function dPdt = modelPertumbuhan(t, P, params)
% MODELPERTUMBUHAN Laju pertumbuhan biomassa model logistik.
%   dPdt = modelPertumbuhan(t, P, params) menghitung turunan dP/dt
%   untuk persamaan logistik:
%       dP/dt = r * P * (1 - P/K)
%   dengan params struct berisi field:
%       r - laju pertumbuhan intrinsik (per hari)
%       K - kapasitas lingkungan (satuan biomassa)
%
%   Argumen t tidak muncul di rumus (model otonom) tetapi wajib
%   diterima karena ode45 selalu memanggil fungsi model dengan
%   pasangan (t, P).
r = params.r;
K = params.K;
dPdt = r * P * (1 - P / K);
end

Mari kita bedah baris-baris kode di atas:

1. Baris 1 (function dPdt = modelPertumbuhan(t, P, params)): tanda tangan fungsi meniru kontrak solver — waktu dulu, keadaan kemudian — lalu menambah argumen ketiga params untuk parameter model.
2. Baris 10-12 (komentar tentang argumen t): t tidak dipakai rumus karena model ini otonom. Tetap menerima t membuat fungsi kompatibel dengan solver tanpa pembungkus tambahan, dan dengan solver non-otonom di proyek lain (misalnya laju musiman di Bagian 4).
3. Baris 13-14 (r = params.r; K = params.K;): parameter diambil dari struct sekali di awal. Seluruh proyek memakai struct params — bukan variabel global — sehingga tiap simulasi membawa fotokopi parameternya sendiri.
4. Baris 15 (dPdt = r * P * (1 - P / K);): seluruh model logistik hidup dalam satu baris. Operasinya skalar demi skalar karena keadaan P di sini memang skalar (satu jenis biomassa).

Bentuk alternatif yang juga sah adalah anonymous function @(t, P) r * P * (1 - P/K). Kita tetap memilih berkas fungsi karena tiga alasan: parameter masuk sebagai argumen sehingga fungsi bisa dipakai ulang oleh sweep tanpa menangkap ulang variabel, fungsi bisa diuji langsung dari tests/, dan blok komentar dokumentasinya menjadi kontrak yang terbaca lewat help modelPertumbuhan.

Menjalankan ode45: jalankanSimulasi.m #

MATLAB menyediakan banyak solver ODE; ode45 adalah kerja kuda untuk model non-stiff dan halus seperti logistik. Solver ini menjalankan Runge-Kutta adaptif orde 4(5) — langkahnya membesar saat kurva landai dan mengecil saat kurva curam — dan sudah kita pakai pada chapter Simulasi Sistem dan Pemodelan Fisika. Berikut seluruh isi src/jalankanSimulasi.m:

function hasil = jalankanSimulasi(params, tspan)
% JALANKANSIMULASI Menjalankan simulasi logistik dengan ode45.
%   hasil = jalankanSimulasi() menjalankan parameter default:
%       r  = 0.08 per hari, K = 100 satuan, P0 = 1,
%       output setiap hari selama 120 hari (tspan = 0:1:120).
%   hasil = jalankanSimulasi(params) memakai params kustom (struct
%       dengan field r, K, P0).
%   hasil = jalankanSimulasi(params, tspan) memakai vektor waktu
%       output sendiri (harus menaik, dimulai dari 0).
%
%   hasil : struct berisi field
%       t              - vektor waktu (hari), 121x1 untuk default
%       P              - biomassa (satuan), sepanjang t
%       params         - parameter yang dipakai simulasi
%       waktuSetengahK - t saat P = K/2 (interpolasi linear, hari);
%                        NaN bila P tidak pernah mencapai K/2 pada
%                        rentang tspan
%
%   Dua keputusan membuat hasil deterministik lintas MATLAB/Octave:
%   tspan diberikan sebagai vektor waktu penuh sehingga kedua mesin
%   wajib mengeluarkan nilai tepat pada waktu itu, dan toleransi
%   diikat ketat lewat odeset. Dua pemanggilan menghasilkan matriks
%   yang identik (isequal).
if nargin < 1 || isempty(params)
    params = struct('r', 0.08, 'K', 100, 'P0', 1);
end
if nargin < 2 || isempty(tspan)
    tspan = 0:1:120;   % output tepat pada hari ke-0 sampai ke-120
end

opsi = odeset('RelTol', 1e-10, 'AbsTol', 1e-12);
[t, P] = ode45(@(t, P) modelPertumbuhan(t, P, params), tspan, ...
    params.P0, opsi);
t = t(:);
P = P(:);

% P monoton naik sehingga interp1(P, t, ...) sah dipakai.
waktuSetengahK = interp1(P, t, params.K / 2);

hasil = struct('t', t, 'P', P, 'params', params, ...
    'waktuSetengahK', waktuSetengahK);
end

Mari kita bedah baris-baris kode di atas:

1. Baris 24-26 (if nargin < 1 ... end): pola optional arguments klasik. jalankanSimulasi() mendapat parameter default; jalankanSimulasi(params) membawa parameternya sendiri — bentuk yang dipakai sweepParameter nanti.
2. Baris 27-29 (tspan = 0:1:120;): inilah keputusan determinisme pertama. tspan dengan lebih dari dua elemen memaksa solver mengeluarkan nilai tepat pada tiap waktu yang diminta — di sini setiap hari ke-0 sampai ke-120, 121 titik — bukan pada waktu langkah internalnya yang bebas dipilih.
3. Baris 31 (opsi = odeset('RelTol', 1e-10, 'AbsTol', 1e-12);): keputusan determinisme kedua. Toleransi default ode45 longgar; mengikatnya ketat membuat tiap langkah internal begitu akurat sehingga nilai yang diinterpolasi ke grid harian konvergen ke solusi eksak sampai jauh di bawah 1e-6 — di kedua mesin.
4. Baris 32-33 (ode45(@(t, P) modelPertumbuhan(t, P, params), ...)): anonymous function membungkus modelPertumbuhan menjadi fungsi dua argumen (t, P) seperti yang diminta solver, sambil menangkap params dari ruang lingkup pemanggil. Kondisi awal diambil dari params.P0.
5. Baris 34-35 (t = t(:); P = P(:);): menegaskan kedua keluaran berbentuk vektor kolom — bentuk yang diandalkan analisisFase dan uji otomatis.
6. Baris 38 (waktuSetengahK = interp1(P, t, params.K / 2);): P monoton naik, jadi interp1 boleh dipakai terbalik — mencari waktu pada biomassa tertentu. Hasilnya t50.
7. Baris 40-41 (hasil = struct(...)): seluruh keluaran dirangkai satu struct — data, parameter, dan metrik — sehingga fungsi-fungsi berikutnya cukup menerima satu objek.

Menjalankannya dari Command Window pada folder src/:

>> hasil = jalankanSimulasi();
>> fprintf('P(120) = %.4f\n', hasil.P(end));
P(120) = 99.3340
>> fprintf('t50 = %.4f hari\n', hasil.waktuSetengahK);
t50 = 57.4390 hari

Validasi terhadap Solusi Analitik #

Kekuatan proyek ini: kita punya jawaban eksak dari Bagian 1, jadi hasil numerik bisa diadili. Potongan berikut diambil langsung dari tests/uji_simulasi.m:

% Solusi analitik P(t) = K/(1 + (K/P0 - 1)*exp(-r*t)) sebagai pembanding.
Panalitik = hasil.params.K ./ ...
    (1 + (hasil.params.K / hasil.params.P0 - 1) * ...
    exp(-hasil.params.r * hasil.t));
galatMaks = max(abs(hasil.P - Panalitik));

Mari kita bedah baris-baris kode di atas:

1. Baris 2-4 (Panalitik = ...): rumus analitik Bagian 1 dievaluasi pada seluruh grid waktu hasil.t sekaligus — operasi vektor elemen demi elemen, tanpa loop.
2. Baris 5 (galatMaks = max(abs(hasil.P - Panalitik));): galat mutlak terbesar di seluruh 121 titik. Satu angka yang merangkum seberapa jujur solver pada model ini.

Uji otomatis kemudian menuntut galatMaks < 1e-6; kasusnya lolos dengan pesan LULUS: galat terhadap solusi analitik < 1e-6. Bagian ringkasan bernilai dari tests/uji_simulasi.m mencetak:

Output yang Dihasilkan:

Ringkasan simulasi default (r = 0.08/hari, K = 100, P0 = 1):
  titik output  : 121 (hari 0 sampai 120)
  P(0)          : 1.0000
  P(120)        : 99.3340
  t50           : 57.4390 hari
  laju maksimum : 1.9997 per hari
  galat maksimum terhadap analitik: 0.000000

Struktur Hasil dan Makna t50 #

Struct hasil membawa empat field yang menjadi bahasa bersama seluruh proyek:

Nilai t50 hasil interp1 adalah 57.4390 hari — berimpit dengan rumus analitik ln(99)/0.08 = 57.4390 sampai empat desimal. Selisih keduanya di bawah 0.001 hari, sisa galat interpolasi linear pada grid harian. Metrik kedua yang dicetak uji, laju maksimum = 1.9997 per hari, hanya 0.0003 di bawah nilai teoretis r*K/4 = 2 — dan itu puncak yang dicari diff(P)./diff(t) pada grid satu hari, bukan turunan eksak.

Sensitivitas Toleransi #

Seberapa penting odeset tadi? Tabel berikut menjalankan simulasi yang sama dengan tiga pengaturan toleransi, dicatat dari Octave:

Dua baris pertama memakai toleransi longgar, dan dua angka pada baris itu (galat dan t50) adalah hasil urutan langkah adaptif solver: MATLAB bisa memilih langkah yang berbeda sehingga digit-digitnya ikut berbeda — orde besarannya tetap sama. Baris ketiga bebas dari masalah itu: toleransi yang diikat membuat hasil konvergen penuh ke analitik, sama di mesin mana pun.

Toleransi default sebenarnya cukup untuk menggambar kurva — galat 6e-3 tak terlihat mata. Tapi proyek ini mencetak angka dan menguji klaim seperti “laju maksimum di P = K/2”, jadi kita butuh angka yang stabil hingga empat desimal; baris ketiga itulah wilayah kerja simulasi ini.

Mesin dan solver akhirnya bekerja dan hasilnya lulus audit analitik. Angka yang bisa dipercaya sudah ada di tangan; Bagian 3 mengubahnya menjadi tiga gambar yang bisa dibaca dalam hitungan detik — kurva, laju, dan perbandingan skenario.