DifferentialEquations.jl

DifferentialEquations.jl adalah meta-package suite SciML untuk pemecahan persamaan diferensial. Topik ini membahas penyelesaian persamaan diferensial biasa (ordinary differential equation, ODE) mulai dari bentuk masalah, konstruksi ODEProblem, pemanggilan solve dengan solver eksplisit Tsit5, pengaturan jendela simpan melalui saveat, hingga event dan callback melalui ContinuousCallback. Seluruh contoh dijalankan dengan OrdinaryDiffEq, yaitu package solver ODE inti yang menjadi anggota utama meta-package tersebut.

Bentuk masalah dan ODEProblem#

Setiap package solver ODE pada suite SciML bekerja pada sistem persamaan diferensial orde satu yang dinyatakan dalam bentuk baku:

Simbol $u$ menyatakan keadaan (state) sistem yang berupa satu bilangan atau array, $p$ menyatakan parameter yang diteruskan ke fungsi tanpa diubah, dan $t$ menyatakan waktu. Fungsi $f$ menghitung laju perubahan keadaan, sedangkan $u_0$ adalah kondisi awal pada waktu $t_0$. Bentuk orde satu tidak membatasi masalah yang dapat diselesaikan karena persamaan orde tinggi selalu dapat dipecah menjadi beberapa persamaan orde satu.

Fungsi dinamika ditulis dengan dua gaya. Gaya in-place, yaitu f(du, u, p, t), menerima array kerja du dan menuliskan hasil ke dalamnya tanpa membuat array baru. Gaya out-of-place, yaitu f(u, p, t), mengembalikan laju perubahan sebagai nilai balik. Konvensi penamaannya menandai perbedaan itu: fungsi in-place diakhiri tanda seru mengikuti kebiasaan Julia untuk fungsi yang mengubah argumennya. Contoh berikut menghitung laju perubahan pertumbuhan logistik secara out-of-place:

# Bentuk out-of-place: f(u, p, t) mengembalikan laju perubahan
dinamika(u, p, t) = p[1] * u * (1 - u / p[2])

println(dinamika(0.5, (0.8, 10.0), 0.0))

0.38

Pertumbuhan logistik digunakan sebagai contoh sepanjang topik karena rumusnya sederhana dan perilakunya mudah dibayangkan:

Parameter $r$ adalah laju pertumbuhan, sedangkan $K$ adalah kapasitas daya dukung yang menahan pertumbuhan ketika populasi mendekatinya. Nilai laju pada keadaan awal $u = 0.5$ adalah 0,38 sebagaimana ditunjukkan keluaran di atas. Masalah lengkap dirangkai dengan ODEProblem(f, u0, tspan, p), yaitu objek yang menyimpan fungsi dinamika, kondisi awal, selang waktu, dan parameter dalam satu wadah.

Menyelesaikan masalah dan membaca solusi#

Pemecahan masalah dilakukan dengan solve(masalah, Tsit5()), yaitu pemanggilan yang menerima objek masalah dan solver yang dipilih. Solver Tsit5 adalah metode Runge-Kutta eksplisit orde lima yang menyesuaikan panjang langkah secara adaptif: langkah diperbesar ketika solusi berubah perlahan dan diperkecil ketika solusi berubah cepat. Hasilnya berupa objek solusi yang menyimpan seluruh lintasan, dengan sol.t berisi waktu tiap titik dan sol.u berisi keadaan pada waktu tersebut.

using OrdinaryDiffEq

# Pertumbuhan logistik: du/dt = r u (1 - u/K)
function logistik!(du, u, p, t)
    r, K = p
    du[1] = r * u[1] * (1 - u[1] / K)
end

u0 = [0.5]
tspan = (0.0, 10.0)
p = (0.8, 10.0)

masalah = ODEProblem(logistik!, u0, tspan, p)
solusi = solve(masalah, Tsit5())

println("Nilai akhir u: ", solusi.u[end][1])
println("Waktu akhir  : ", solusi.t[end])
println("Jumlah titik tersimpan: ", length(solusi.t))

Nilai akhir u: 9.936579267188284
Waktu akhir  : 10.0
Jumlah titik tersimpan: 13

Nilai akhir mendekati kapasitas daya dukung 10.0 sesuai perilaku logistik, dan selisih kecil pada angka keempat masih berada dalam rentang toleransi bawaan solver. Jumlah titik tersimpan tidak ditentukan sebelumnya karena mengikuti langkah adaptif. Keadaan pada titik mana pun diambil lewat indeks, misalnya sol.u[5] untuk titik kelima, dan objek solusi juga dapat dipanggil dengan argumen waktu sehingga sol(2.5) mengembalikan interpolasi keadaan pada waktu tersebut. Alur dari fungsi dinamika hingga objek solusi dirangkum pada diagram berikut.

Jendela simpan (save window) mengatur titik waktu mana yang disimpan ke dalam objek solusi. Argumen saveat menerima vektor waktu, dan nilai pada titik yang bukan ujung langkah diperoleh dari interpolasi. Penyimpanan yang teratur memudahkan penyusunan tabel dan grafik:

solusi_tersampel = solve(masalah, Tsit5(); saveat = 0.0:2.5:10.0)
println(solusi_tersampel.t)
for (t, u) in zip(solusi_tersampel.t, solusi_tersampel.u)
    println("t = ", t, "  ->  u = ", round(u[1], digits = 3))
end

[0.0, 2.5, 5.0, 7.5, 10.0]
t = 0.0  ->  u = 0.5
t = 2.5  ->  u = 2.8
t = 5.0  ->  u = 7.418
t = 7.5  ->  u = 9.55
t = 10.0  ->  u = 9.937

Event dan callback dengan ContinuousCallback#

Event adalah peristiwa pada waktu yang tidak diketahui sebelumnya, misalnya saat populasi menyentuh sebuah ambang. Suite SciML menangani event melalui callback, yaitu fungsi yang dipanggil solver pada titik tertentu. ContinuousCallback memantau sebuah fungsi kondisi dan mendeteksi perlintasan nolnya dengan pencarian akar, lalu menjalankan fungsi affect saat kejadian itu terjadi. Apabila fungsi affect memanggil terminate!, integrasi berhenti tepat di titik kejadian.

Mekanisme deteksinya bekerja di antara dua langkah. Solver mengamati tanda fungsi kondisi pada kedua ujung langkah, dan perubahan tanda menandakan perlintasan nol yang kemudian dicari letak persisnya. Cara ini membuat waktu kejadian tidak terikat pada panjang langkah, sehingga hasilnya presisi meskipun langkah adaptifnya besar. Konstruktor yang sama juga menerima fungsi affect kedua untuk arah perlintasan sebaliknya, yaitu ketika kondisi turun melewati nol, sehingga naik dan turun dapat ditangani secara berbeda.

using OrdinaryDiffEq

dinamika(u, p, t) = p[1] * u * (1 - u / p[2])

masalah = ODEProblem(dinamika, 0.5, (0.0, 20.0), (0.8, 10.0))

# Kondisi bernilai nol tepat saat populasi mencapai 8.0
kondisi(u, t, integrator) = u - 8.0
cb = ContinuousCallback(kondisi, integrator -> terminate!(integrator))

solusi = solve(masalah, Tsit5(); callback = cb)

println("Populasi saat integrasi berhenti: ", round(solusi.u[end], digits = 3))
println("Waktu kejadian: ", round(solusi.t[end], digits = 4))

Populasi saat integrasi berhenti: 8.0
Waktu kejadian: 5.4134

Integrasi berhenti pada waktu 5,4134 dengan populasi tepat 8,0, padahal tspan mencakup hingga waktu 20.0. Callback juga dapat dipakai untuk pola lain: fungsi affect boleh mengubah keadaan melalui integrator.u, misalnya membalik arah pertumbuhan, tanpa menghentikan integrasi. Satu masalah dapat membawa beberapa callback sekaligus melalui CallbackSet, sehingga kejadian yang berbeda jenisnya ditangani dalam satu pemanggilan solve.

Berkas contoh: pendinginan.jl#

Pendinginan Newton memodelkan benda yang suhunya $T$ mendekati suhu lingkungan:

Konstanta $k$ mengatur cepatnya pendinginan, dan kasus berikut mengikuti catatan sensor suhu dengan suhu awal 100 °C serta suhu lingkungan 25 °C. Berkas pendinginan.jl merupakan program utuh yang merumuskan masalah tersebut, menyelesaikannya dengan Tsit5, lalu menampilkan suhu akhir beserta jumlah titik waktu tersimpan:

# pendinginan.jl
# Menyelesaikan pendinginan Newton dengan OrdinaryDiffEq.

using OrdinaryDiffEq

# Pendinginan Newton: dT/dt = -k (T - T_lingkungan)
function pendinginan!(dT, T, p, t)
    k, T_lingkungan = p
    dT[1] = -k * (T[1] - T_lingkungan)
end

k = 0.07
T_lingkungan = 25.0
T_awal = [100.0]

masalah = ODEProblem(pendinginan!, T_awal, (0.0, 60.0), (k, T_lingkungan))
solusi = solve(masalah, Tsit5())

println("Suhu pada t = ", solusi.t[end], " s: ",
        round(solusi.u[end][1], digits = 2), " °C")
println("Jumlah titik waktu tersimpan: ", length(solusi.t))

Suhu pada t = 60.0 s: 26.12 °C
Jumlah titik waktu tersimpan: 11

Suhu akhir 26,12 °C konsisten dengan solusi analitik $25 + 75e^{-4.2} \approx 26{,}12$. Konsistensi seperti ini dimungkinkan oleh sifat persamaannya yang linear, sehingga masalah pendinginan Newton sering menjadi uji kecil sebelum solver diterapkan pada model yang rumit. Selisih suhu terhadap lingkungan menurun secara eksponensial, dan dengan $k = 0.07$ selisih 75 °C pada awal pengamatan tinggal sekitar 1,1 °C setelah 60 detik. Berkas semacam ini berfungsi sebagai template karena bagian yang perlu diubah untuk masalah lain hanya fungsi dinamika, kondisi awal, dan parameter. Programnya dijalankan dari direktori topik dengan perintah julia pendinginan.jl.

Kesalahan umum#

Kesalahan yang paling sering muncul saat mempraktikkan topik ini menyangkut fungsi dinamika dan callback. Daftar berikut memuat penyebab dan perbaikannya.

• Tidak mengisi seluruh elemen du pada fungsi in-place. Array du diterima sebagai memori kerja yang belum diinisialisasi, sehingga elemen yang tidak ditulis berisi nilai acak dan hasil integrasi salah tanpa pesan error apa pun. Perbaikannya dengan menulis setiap elemen du atau beralih ke bentuk out-of-place.
• Menganggap callback otomatis menghentikan integrasi. Fungsi affect hanya menjalankan yang diperintahkan, sehingga integrasi berlanjut sampai ujung tspan apabila terminate! tidak dipanggil.
• Menyatakan panjang keadaan yang tidak cocok dengan dimensi operasi di dalam fungsi dinamika. Error yang muncul biasanya DimensionMismatch, dan contoh berikut menangkapnya dengan try-catch agar program selesai dengan normal:

using OrdinaryDiffEq, LinearAlgebra

A = [1.0 2.0; 3.0 4.0]
u0 = [1.0]  # panjang keadaan tidak cocok dengan dimensi A

f_salah!(du, u, p, t) = mul!(du, A, u)

try
    solve(ODEProblem(f_salah!, u0, (0.0, 1.0)), Tsit5())
catch e
    println("Tertangkap: ", typeof(e))
end

Tertangkap: DimensionMismatch

Pesan DimensionMismatch muncul saat solver memanggil f_salah! pertama kali karena mul! menolak perkalian matriks 2×2 dengan vektor yang hanya berisi satu elemen. Membaca tipe errornya sudah cukup untuk menemukan penyebabnya, yaitu panjang u0 yang tidak sesuai dengan dimensi kolom matriks A pada perkalian A * u.

Ringkasan#

• Masalah ODE dinyatakan dengan ODEProblem(f, u0, tspan, p) dalam bentuk orde satu, dengan fungsi dinamika in-place f(du, u, p, t) atau out-of-place f(u, p, t).
• solve(masalah, Tsit5()) menghasilkan objek solusi dengan sol.t dan sol.u, sedangkan saveat mengatur titik waktu yang disimpan melalui interpolasi.
• ContinuousCallback mendeteksi perlintasan nol fungsi kondisi lalu menjalankan fungsi affect, dan pemanggilan terminate! menghentikan integrasi tepat pada titik kejadian.

Referensi#

• DifferentialEquations.jl Documentation — SciML — Dokumentasi resmi meta-package pemecah persamaan diferensial beserta tutorial dan pilihan solver.
• SciMLBase Documentation — SciML — Definisi tipe masalah seperti ODEProblem dan antarmuka solve yang dipakai seluruh package SciML.