Pemodelan Fisika

Fisika memberi pemodelan komputasi bonus yang jarang ada di bidang lain: hukum yang sudah teruji ratusan tahun dan solusi analitik yang bisa jadi hakim. Ketika simulasi numerik kita menyimpang dari rumus parabola atau solusi osilasi teredam, yang salah hampir pasti kode atau modelnya — bukan alam semesta. Chapter ini menggarap tiga kasus klasik dengan tiga jalur penyelesaian berbeda: gerak parabola sepenuhnya analitik, osilasi teredam via ode45 yang diadukan ke rumusnya, dan pendinginan Newton sebagai model kalor, semuanya ditutup verifikasi kekekalan energi.

• Gerak parabola: solusi analitik dan jarak maksimum
• Osilasi teredam m x'' + c x' + k x = 0 dengan vektor state [x; v]
• Verifikasi energi: awal = sisa + terdisipasi (integral trapz)
• Pendinginan Newton dT/dt = -k(T - Tenv)

Dari Hukum Fisika ke Persamaan #

Alur kerjanya selalu sama: tulis hukumnya, pilih variabel, selesaikan. Yang sering menjegal bukan matematikanya, melainkan satuan: persamaan hanya konsisten bila seluruh besaran dibawa ke satu sistem satuan (kita pakai SI — meter, sekon, kilogram, joule). Tiga kasus chapter ini memetakan tiga strategi penyelesaian:

Gerak Parabola: Solusi Analitik #

Proyektil yang diluncurkan dengan kecepatan v0 dan elevasi sudut membelah menjadi dua komponen: horizontal v0x (konstan, tanpa gaya) dan vertikal v0y (melawan gravitasi). Seluruh pertanyaan penting — waktu terbang, jarak, tinggi maksimum — punya rumus tertutup:

g     = 9.81;              % gravitasi (m/s^2)
v0    = 20;                % kecepatan awal (m/s)
sudut = 35*pi/180;         % elevasi (radian)
v0x = v0*cos(sudut);       % komponen horizontal
v0y = v0*sin(sudut);       % komponen vertikal
tTerbang = 2*v0y/g;
jarak    = v0x*tTerbang;
tinggi   = v0y^2/(2*g);
fprintf('parabola: t terbang = %.4f s, jarak = %.4f m, tinggi maks = %.4f m\n', ...
        tTerbang, jarak, tinggi);
fprintf('jarak rumus v0^2*sin(2*sudut)/g = %.4f m (identik)\n', ...
        v0^2*sin(2*sudut)/g);
tt = linspace(0, tTerbang, 50);
xx = v0x*tt;
yy = v0y*tt - 0.5*g*tt.^2;
[yPuncak, iPuncak] = max(yy);
fprintf('puncak grid: %.2f m pada t = %.2f s\n', yPuncak, tt(iPuncak));

Output yang Dihasilkan:

parabola: t terbang = 2.3387 s, jarak = 38.3157 m, tinggi maks = 6.7072 m
jarak rumus v0^2*sin(2*sudut)/g = 38.3157 m (identik)
puncak grid: 6.70 m pada t = 1.15 s

Mari kita bedah baris-baris kode di atas:

1. Baris 3 (35*pi/180): trigonometri MATLAB memakai radian — elevasi 35 derajat harus dikonversi lebih dulu.
2. Baris 6-8 (tTerbang, jarak, tinggi): waktu terbang dari syarat y = 0; jarak = v0x kali waktu terbang; tinggi dari kekekalan energi kinetik-vertikal menjadi potensial.
3. Baris 13-16 (linspace … max): 50 titik grid waktu menghasilkan lintasan mulus; max mengembalikan nilai puncak beserta indeksnya.
4. Baris 7 dan 11-12 (verifikasi silang): jarak via v0x*tTerbang dan via v0^2*sin(2*sudut)/g menghasilkan 38.3157 m — dua rumus berbeda, satu angka.

Lintasan xx-yy yang di-plot dengan axis equal menggambar parabola simetris; contoh.m menyimpan grafiknya sebagai PNG lewat pola headless biasa.

Osilasi Teredam: ode45 vs Solusi Analitik #

Massa-pegas-damper adalah sistem uji paling jujur untuk solver ODE. Persamaan m*x'' + c*x' + k*x = 0 dipecah menjadi pasangan orde satu dengan state [x; v], lalu solusi numeriknya dibandingkan dengan rumus underdamped. Mekanisme pembandingnya: rasio redaman zeta = c/(2*sqrt(k*m)) memutuskan watak sistem — zeta < 1 berarti osilasi meluruh, dan solusi analitiknya adalah exp(-zeta*w0*t) dikali cosinus pada frekuensi teredam wd:

m = 1; c = 0.8; k = 16;    % massa, redaman, pegas (SI)
w0   = sqrt(k/m);
zeta = c/(2*sqrt(k*m));
wd   = w0*sqrt(1 - zeta^2);
fprintf('osilasi: w0 = %.4f rad/s, zeta = %.4f, wd = %.4f rad/s\n', w0, zeta, wd);
osilasi = @(t,y) [y(2); (-c*y(2) - k*y(1))/m];
opt = odeset('RelTol', 1e-10, 'AbsTol', 1e-12);
[t, y] = ode45(osilasi, [0 1 2 3 4 5], [0.05; 0], opt);
xAn = @(t) 0.05*exp(-zeta*w0*t) .* (cos(wd*t) + (zeta*w0/wd)*sin(wd*t));
fprintf('%4s %12s %12s\n', 't', 'x_ode45', 'x_analitik');
for i = 1:numel(t)
    fprintf('%4d %12.6f %12.6f\n', t(i), y(i,1), xAn(t(i)));
end

Output yang Dihasilkan:

osilasi: w0 = 4.0000 rad/s, zeta = 0.1000, wd = 3.9799 rad/s
   t      x_ode45   x_analitik
   0     0.050000     0.050000
   1    -0.024916    -0.024916
   2    -0.000130    -0.000130
   3     0.011312     0.011312
   4    -0.010082    -0.010082
   5     0.003956     0.003956

Mari kita bedah baris-baris kode di atas:

1. Baris 2-4 (w0, zeta, wd): parameter fisika diringkas dulu — frekuensi natural, rasio redaman, frekuensi teredam — supaya rumus analitik terbaca sebagai rumus bukan sebagai tumpukan konstanta.
2. Baris 6 (osilasi = @(t,y) ...): persamaan orde dua ditulis ulang x' = v, v' = (-c*v - k*x)/m — persis pola vektor state chapter simulasi-sistem.
3. Baris 8 ([0 1 2 3 4 5]): rentang enam elemen meminta output pada waktu-waktu kunci; toleransi 1e-10 memastikan digit ke-enam bisa dipercaya.
4. Baris 9 (xAn): rumus underdamped untuk kondisi awal x0 = 0.05, v0 = 0 — suku sin ekstra muncul dari kecepatan awal nol.

Kolom x_ode45 dan x_analitik identik sampai enam desimal pada keenam titik waktu. Simulasi digital dan rumus 1850-an sepakat.

Verifikasi Kekekalan Energi #

Sistem teredam melanggar kekekalan energi mekanik — damper menyerapnya. Hukum yang tetap berlaku: energi awal = energi tersisa + energi terdisipasi, dan daya disipasi damper adalah c*v^2. Integralnya dihitung trapz dari lintasan kecepatan:

[tD, yD] = ode45(osilasi, [0 8], [0.05; 0], opt);
xD = yD(:,1);  vD = yD(:,2);
E    = 0.5*m*vD.^2 + 0.5*k*xD.^2;      % energi mekanik sesaat
Edis = c*trapz(tD, vD.^2);             % energi serap damper
fprintf('energi awal   = %.6f J\n', E(1));
fprintf('energi sisa   = %.6f J\n', E(end));
fprintf('energi disipasi = %.6f J\n', Edis);
selisih = abs(E(1) - (E(end) + Edis))/E(1);
if selisih < 1e-5
    fprintf('keseimbangan energi: LULUS (selisih relatif < 1e-5)\n');
else
    fprintf('keseimbangan energi: GAGAL (selisih relatif %.2e)\n', selisih);
end

Output yang Dihasilkan:

energi awal   = 0.020000 J
energi sisa   = 0.000036 J
energi disipasi = 0.019964 J
keseimbangan energi: LULUS (selisih relatif < 1e-5)

Energi awal 0.02 J terurai rapi: 0.000036 J masih tersimpan sebagai osilasi sisa, 0.019964 J terserap damper — buku kasnya seimbang hingga galat relatif di bawah satu per seratus ribu. Verifikasi semacam ini tidak butuh solusi analitik; ia hanya butuh hukum kekekalan, sehingga tetap berlaku untuk sistem yang jauh lebih rumit.

Pendinginan Newton #

Kopi 90 derajat di ruang 25 derajat mendingin dengan laju sebanding selisih suhunya: dT/dt = -k*(T - Tenv). Model orde satu ini adalah kasus uji ode45 paling sederhana sekaligus contoh nyata peluruhan eksponensial menuju kesetimbangan:

kN = 0.08; Tenv = 25; T0 = 90;   % konstanta, suhu ruang, suhu awal
pendingin = @(t,T) -kN*(T - Tenv);
[tT, TT] = ode45(pendingin, [0 10 20 30 45], T0, opt);
TAn = @(t) Tenv + (T0 - Tenv)*exp(-kN*t);
fprintf('%4s %12s %12s\n', 't', 'T ode45', 'T analitik');
for i = 1:numel(tT)
    fprintf('%4d %12.4f %12.4f\n', tT(i), TT(i), TAn(tT(i)));
end
fprintf('waktu mencapai 40 C = %.4f s\n', log((T0-Tenv)/(40-Tenv))/kN);

Output yang Dihasilkan:

   t      T ode45   T analitik
   0      90.0000      90.0000
  10      54.2064      54.2064
  20      38.1233      38.1233
  30      30.8967      30.8967
  45      26.7760      26.7760
waktu mencapai 40 C = 18.3292 s

Baris terakhir membalik modelnya: rumus eksponensial diselesaikan untuk T = 40 dengan sekali log, menghasilkan 18.33 sekon tanpa perlu mencari akar secara numerik. Suhu menuju Tenv = 25 namun tak pernah tiba — pada 45 sekon masih tersisa selisih 1.78 derajat.

Berkas contoh.m menjalankan ketiga model berurutan, memverifikasi kesetimbangan energi, dan menyimpan plot lintasan parabola — dalam satu script yang selesai dalam hitungan detik.

Ringkasan #

Tiga kasus, tiga strategi: parabola diselesaikan rumus tertutup (jarak = 38.3157 m, terkonfirmasi dua rumus berbeda), osilasi teredam diselesaikan ode45 dengan vektor state [x; v] dan cocok enam desimal dengan rumus underdamped, pendinginan Newton dijalankan ganda — numerik dan analitik — plus inversi logaritma untuk waktu pendinginan 18.3292 s. Verifikasi kekekalan energi menjembatani keduanya: 0.02 J awal terbagi menjadi 0.000036 J sisa plus 0.019964 J disipasi via trapz, dengan selisih relatif di bawah 1e-5. Satu disiplin berdiri di seluruhnya: satuan SI yang konsisten, toleransi ketat untuk angka yang dibandingkan, dan hukum kekekalan sebagai hakim yang tidak pernah menyogok.