Simulasi Sistem

Kadang kita tidak butuh meng-fit data, melainkan menjawab pertanyaan “jika”: jika populasi dibiarkan 20 jam, berapa jumlahnya? jika pegas dilepas dari posisi 10 cm, di mana massa itu tiga detik kemudian? Simulasi adalah eksperimen virtual — kita menulis hukum perubahan sistem sebagai persamaan diferensial, lalu membiarkan komputer mengintegrasikannya langkah demi langkah. MATLAB membungkus mekanisme numeriknya dalam satu panggilan: ode45. Chapter ini menyusun model dengan fungsi anonim, menjalankannya untuk dua sistem (pertumbuhan populasi dan osilasi pegas), lalu memvalidasi hasilnya lewat hukum kekekalan energi.

• Persamaan diferensial orde satu dan fungsi anonim sebagai model
• Panggilan dasar ode45: rentang waktu, kondisi awal, hasil [t, y]
• Sistem dua persamaan dengan vektor state [x; v]
• Toleransi via odeset dan validasi kekekalan energi

Sistem Dinamik sebagai Persamaan Diferensial #

Aturan simulasi dinamik: nyatakan laju perubahan, bukan nilainya. Sistem dinamik ditulis sebagai dy/dt = f(t, y) — laju perubahan state y ditentukan oleh waktu dan state itu sendiri. Populasi yang tumbuh makin lambat seiring menuju kapasitas lingkungan, misalnya, mengikuti model logistik dP/dt = r*P*(1 - P/K): saat P masih kecil, lajunya hampir r*P; saat P mendekati K, laju menuju nol. Mekanisme di balik layar: solver tidak bisa melompat ke solusi langsung — solver menghitung f(t, y) di titik sekarang, maju sedikit, hitung ulang, dan seterusnya, sehingga hasil akhirnya adalah tabel (t, y) hasil integrasi numerik bertahap.

ode45: Menjalankan Model dari Fungsi Anonim #

Model kita hidup sebagai fungsi anonim — persamaan diferensial dibungkus handle @(t, y) yang mengembalikan turunan. Fungsi ode45 menerima handle itu, rentang waktu, dan kondisi awal:

r  = 0.35;    % laju pertumbuhan intrinsik (per jam)
K = 500;      % kapasitas lingkungan (individu)
P0 = 20;      % populasi awal
logistik = @(t,y) r*y*(1 - y/K);
opt = odeset('RelTol', 1e-10, 'AbsTol', 1e-12);

[tLintasan, P]  = ode45(logistik, [0 20], P0, opt);       % lintasan padat
[tKunci, PK]    = ode45(logistik, [0 5 10 20], P0, opt);  % output waktu kunci
analitik = @(t) K ./ (1 + (K/P0 - 1)*exp(-r*t));
fprintf('%4s %10s %10s\n', 't', 'P ode45', 'P analitik');
for i = 1:numel(tKunci)
    fprintf('%4d %10.4f %10.4f\n', tKunci(i), PK(i), analitik(tKunci(i)));
end

Output yang Dihasilkan:

   t    P ode45 P analitik
   0    20.0000    20.0000
   5    96.7011    96.7011
  10   289.8992   289.8992
  20   489.2918   489.2918

Mari kita bedah baris-baris kode di atas:

1. Baris 4 (logistik = @(t,y) ...): fungsi anonim menampung seluruh model — argumen pertama selalu waktu t, argumen kedua selalu state y, dan nilai kembalian adalah dy/dt.
2. Baris 5 (odeset): option set mengatur toleransi. RelTol membatasi galat relatif per langkah, AbsTol galat absolut untuk nilai kecil.
3. Baris 7 ([0 20]): rentang waktu dua elemen — solver memilih langkah adaptif sendiri dan mengembalikan lintasan padat untuk plot.
4. Baris 8 ([0 5 10 20]): rentang lebih dari dua elemen — solver tetap berintegrasi adaptif, tetapi output dikembalikan persis pada waktu-waktu yang diminta.
5. Baris 9-13: hasil dibandingkan dengan solusi analitik logistik; keempat waktu kunci cocok sampai empat desimal.

Populasi merambat dari 20 menuju 489.29 — mendekati kapasitas K = 500 tanpa pernah melewatinya, persis seperti yang dijanjikan model.

Sistem Dua Persamaan: Vektor State #

Osilasi pegas melibatkan percepatan: m*x'' = -k*x. Solver ode45 hanya menerima persamaan orde satu, jadi persamaan orde dua dipecah menjadi dua persamaan orde satu dengan state gabungan y = [x; v]:

m  = 1;       % massa (kg)
kp = 16;      % konstanta pegas (N/m)
pegas = @(t,y) [y(2); -kp/m*y(1)];        % y(1)=x, y(2)=v
[tO, yO] = ode45(pegas, [0 10], [0.1; 0], opt);
x = yO(:,1);   v = yO(:,2);
[tK, yK] = ode45(pegas, [0 2.5 5 7.5 10], [0.1; 0], opt);
fprintf('%6s %12s %12s\n', 't', 'x (m)', 'v (m/s)');
for i = 1:numel(tK)
    fprintf('%6.1f %12.6f %12.6f\n', tK(i), yK(i,1), yK(i,2));
end

Output yang Dihasilkan:

     t        x (m)      v (m/s)
   0.0     0.100000     0.000000
   2.5    -0.083907     0.217608
   5.0     0.040808    -0.365178
   7.5     0.015425     0.395213
  10.0    -0.066694    -0.298045

Mari kita bedah baris-baris kode di atas:

1. Baris 3 (pegas = @(t,y) [y(2); ...]): baris pertama keluaran berarti x' = v; baris kedua berarti v' = -k/m*x — Hukum Hooke ditulis ulang sebagai pasangan orde satu.
2. Baris 4 ([0.1; 0]): kondisi awal adalah vektor kolom — posisi 0.1 m, kecepatan nol. Urutan komponen harus konsisten dengan definisi pegas.
3. Baris 5 (yO(:,1), yO(:,2)): hasil yO berukuran banyak-baris × 2-kolom; kolom 1 menyimpan x, kolom 2 menyimpan v — satu kolom per komponen state.
4. Baris 7-10: waktu kunci diminta lewat rentang [0 2.5 5 7.5 10] sehingga tabel tercetak tanpa interpolasi tambahan.

Nilai x berganti tanda secara berkala — osilasi — dan kecepatan v berlawanan fase dengan posisi, persis pola fisika pegas.

Membaca Hasil dan Validasi Energi #

Sistem pegas tanpa redaman menyimpan energi mekanik E = ½mv² + ½kx² yang harus konstan. Solusi analitik tidak selalu tersedia untuk sistem yang kita hadapi, tetapi hukum kekekalan seperti ini selalu bisa dijadikan alat uji:

E = 0.5*m*v.^2 + 0.5*kp*x.^2;
drift = max(abs(E - E(1)))/E(1);
fprintf('energi awal = %.6f J\n', E(1));
if drift < 1e-6
    fprintf('validasi energi: LULUS (drift relatif < 1e-6)\n');
else
    fprintf('validasi energi: GAGAL (drift relatif %.2e)\n', drift);
end

Output yang Dihasilkan:

energi awal = 0.080000 J
validasi energi: LULUS (drift relatif < 1e-6)

Energi awal 0.08 J (seluruhnya energi pegas ½kx²) tidak bocor lebih dari satu per sejuta sepanjang 10 detik simulasi — indikasi bahwa toleransi 1e-10 memadai untuk sistem ini. Validasi model dan validasi numerik berjalan bersamaan: yang pertama menguji apakah persamaan kita benar, yang kedua menguji apakah solver menghormati hukum fisikanya.

Memilih Rentang Waktu dan Toleransi #

Bagaimana jika toleransi default dipakai? Perbandingan berikut menjalankan model pegas yang sama tanpa odeset:

[tDef, yDef] = ode45(pegas, [0 10], [0.1; 0]);            % RelTol default
EDef = 0.5*yDef(:,2).^2 + 0.5*kp*yDef(:,1).^2;
driftDef = max(abs(EDef - EDef(1)))/EDef(1);
fprintf('jumlah langkah : default %d, ketat %d\n', numel(tDef), numel(tO));
fprintf('drift energi   : default %.2e, ketat %.2e\n', driftDef, drift);

Output yang Dihasilkan:

jumlah langkah : default 201, ketat 5845
drift energi   : default 2.37e-03, ketat 4.76e-10

Toleransi default (RelTol = 1e-3) menghasilkan drift energi 0.24 persen — cukup untuk melihat bentuk osilasi, tidak cukup untuk audit energi. Mengetatkan toleransi membayar dengan 29 kali lebih banyak langkah. Pilihan praktisnya: default untuk eksplorasi, ketat untuk hasil final dan validasi hukum kekekalan.

Plot lintasan menutup alur simulasi — pola headless yang sama dengan chapter visualisasi:

fg = figure('visible', 'off');
plot(tO, x, 'b-', 'LineWidth', 1.5);
xlabel('t (s)'); ylabel('x (m)');
title('Osilasi pegas: hasil ode45');
grid on;
saveas(fg, 'lintasan-pegas.png');
fprintf('grafik tersimpan: lintasan-pegas.png\n');
delete(fg);

Grafik x terhadap t menampilkan gelombang sinus amplitudo konstan — bukti visual bahwa sistem tanpa redaman tidak kehilangan energi. Berkas contoh.m menjalankan kedua sistem, validasi energi, perbandingan toleransi, dan plot dalam satu script.

Ringkasan #

Simulasi sistem dimulai dari persamaan dy/dt = f(t, y) yang dibungkus fungsi anonim, lalu ode45 mengintegrasikannya: rentang dua elemen memberi lintasan padat, rentang lebih panjang memberi output persis pada waktu kunci. Sistem orde dua seperti pegas dipecah menjadi vektor state [x; v], dan hasil [t, y] dibaca kolom demi kolom — kolom waktu dulu, lalu satu kolom per komponen state. Validasi tidak menunggu solusi analitik: hukum kekekalan energi memberi uji internal yang murah, dan drift 4.76e-10 pada toleransi 1e-10 membuktikan solver bekerja sesuai fisika. Toleransi default cocok untuk eksplorasi; toleransi ketat adalah harga yang layak dibayar untuk hasil final.