Persamaan diferensial biasa

Menyelesaikan persamaan diferensial biasa dengan pemecah ode45 dan pengaturan toleransi.

Daftar isi

Banyak hukum alam dinyatakan sebagai laju perubahan, dan MATLAB mengerjakan model semacam itu melalui keluarga pemecah persamaan diferensial biasa dengan antarmuka yang seragam. Topik ini membahas bentuk masalah nilai awal, penulisan model sebagai fungsi yang menerima waktu dan keadaan, serta pemecah ode45 beserta struktur keluarannya. Sistem orde dua seperti osilator teredam diturunkan menjadi dua persamaan orde satu, dan toleransi pemecah diatur lewat odeset. Setelah membaca topik ini, model laju perubahan dapat diubah menjadi kurva solusi numerik dalam beberapa baris kode beserta kendali mutunya.

Bentuk masalah nilai awal#

Masalah yang diselesaikan pemecah berbentuk masalah nilai awal,

\frac{dy}{dt} = f(t, y), \qquad y(t_0) = y_0

yaitu laju perubahan $y$ ditentukan oleh fungsi $f$ yang bergantung pada waktu $t$ dan nilai $y$ itu sendiri, kemudian nilai awal $y_0$ menutup satu-satunya kebebasan yang tersisa. Pemecah ode45 mengerjakannya dengan metode Runge-Kutta orde empat-lima yang langkahnya menyesuaikan diri, artinya daerah yang berubah cepat dihitung dengan langkah kecil sedangkan daerah yang datar dihitung dengan langkah besar. Nama ode45 menyimpan kedua orde metodenya, dan pemecah ini menjadi pilihan pertama untuk mayoritas masalah karena keseimbangan kecepatan dan akurasinya.

Keluaran pemecah selalu berpasangan [t, y]. Vektor kolom t memuat titik-titik waktu tempat solusi dievaluasi, dan matriks y memuat satu kolom untuk tiap variabel keadaan (state) sehingga barisnya sejajar dengan t. Banyaknya baris menyatakan banyaknya langkah yang diterima pemecah, dan informasi itu menjadi ukuran biaya komputasi yang deterministik, berbeda dengan waktu eksekusi yang bergantung komputer.

Model sebagai fungsi dan struktur keluarannya#

Model paling ringkas ditulis sebagai fungsi anonim satu baris. Hukum pendinginan Newton, yaitu laju penurunan suhu sebanding dengan selisih suhu terhadap lingkungannya,

\frac{dT}{dt} = -k\,(T - T_{\text{aset}})

diterjemahkan menjadi @(t, T) -k*(T - T_aset) dengan konstanta $k = 0{,}08$ per menit dan suhu lingkungan 20 derajat Celsius. Daftar argumen fungsi model selalu waktu dahulu baru keadaan, walaupun rumusnya tidak menggunakan waktu, karena pemecah memanggil model dengan urutan itu.

matlab Pendinginan Newton dengan fungsi anonim

k = 0.08;                               % per menit
T_aset = 20;                            % derajat Celsius
pendinginan = @(t, T) -k*(T - T_aset);

[t, T] = ode45(pendinginan, [0 60], 90);
fprintf('ukuran t: %s, ukuran T: %s\n', mat2str(size(t)), mat2str(size(T)));
fprintf('banyak langkah: %d\n', numel(t));
fprintf('T pada t=0  : %7.3f\n', T(1));
fprintf('T pada t=60 : %7.3f\n', T(end));

text Keluaran

ukuran t: [81 1], ukuran T: [81 1]
banyak langkah: 81
T pada t=0  :  90.000
T pada t=60 :  20.576

Objek panas 90 derajat mendingin menuju 20 derajat, dan setelah 60 menit suhunya 20.576 derajat, bersisa $70\,e^{-4.8}$ dari selisih awalnya sesuai jawaban analitis $T(t) = 20 + 70\,e^{-0.08t}$ . Kedua keluaran berukuran $(81 \times 1)$ karena satu variabel keadaan menghasilkan satu kolom, dan 81 langkah diterima pemecah pada rentang 60 menit. Penggambaran kurvanya cukup menggunakan plot(t, T) karena t dan T sudah sejajar baris demi baris.

Osilator teredam sebagai sistem orde dua#

Persamaan orde dua, misalnya osilator teredam

y'' + 0.4\,y' + 4\,y = 0

tidak dapat langsung dimasukkan ke pemecah yang hanya mengenal turunan pertama. Cara bakunya adalah menurunkan orde dengan variabel tambahan, yaitu $y_1 = y$ sebagai posisi dan $y_2 = y'$ sebagai kecepatan, sehingga satu persamaan orde dua menjadi dua persamaan orde satu,

y_1' = y_2, \qquad y_2' = -0.4\,y_2 - 4\,y_1

Model sistem seperti ini lebih jelas ditulis sebagai berkas fungsi tersendiri agar badannya dapat memuat beberapa baris. Berkas berikut mendefinisikan model yang menerima t dan vektor keadaan y lalu mengembalikan vektor kolom turunan, dan blok if nargin < 2 memasang masukan contoh agar berkas tetap selesai tanpa error ketika dijalankan tanpa argumen.

matlab osilator_teredam.m

function dydt = osilator_teredam(t, y)
% OSILATOR_TEREDAM model osilator teredam sebagai sistem dua persamaan orde satu.
%   dydt = OSILATOR_TEREDAM(t, y) mengembalikan turunan
%   [y(2); -0.4*y(2) - 4*y(1)] untuk keadaan y = [posisi; kecepatan].
%   Tanpa argumen, contoh t = 0 dan y = [2; 0] dipakai.

if nargin < 2
  t = 0;
  y = [2; 0];
end
dydt = [y(2); -0.4*y(2) - 4*y(1)];
end

Skrip pemanggilnya memecahkan sistem dengan awalan posisi 2 dan kecepatan 0 selama 20 detik, mencetak ringkasannya, lalu menggambar kedua kolom keadaan sebagai kurva. Operator @ pada @osilator_teredam membentuk handle yang menunjuk fungsi pada berkas tersebut, dan pemecah memanggilnya berulang kali di sepanjang rentang waktu.

matlab simulasi_osilator.m

% simulasi_osilator.m: memecah osilator teredam lalu menggambar kurvanya
%
% persamaan: y'' + 0.4*y' + 4*y = 0, dengan y(0) = 2 dan y'(0) = 0

[t, y] = ode45(@osilator_teredam, [0 20], [2; 0]);

fprintf('banyak langkah : %d\n', numel(t));
fprintf('ukuran y       : %s\n', mat2str(size(y)));
fprintf('posisi terakhir: %8.4f\n', y(end, 1));

plot(t, y(:, 1), 'b-', t, y(:, 2), 'r--');
xlabel('waktu (detik)');
ylabel('nilai keadaan');
title('Osilator teredam: posisi dan kecepatan');
legend('posisi y1', 'kecepatan y2', 'location', 'northeast');
print('-dsvg', 'osilator-teredam.svg');

text Keluaran

banyak langkah : 217
ukuran y       : [217 2]
posisi terakhir:  -0.0154

Kurva posisi dan kecepatan hasil skrip itu tampak pada figur berikut.

Solusi osilator teredam: posisi bergaya garis padat dan kecepatan bergaya putus-putus

Matriks y berukuran $(217 \times 2)$ karena dua variabel keadaan menghasilkan dua kolom, dan 217 langkah diterima pemecah, lebih banyak daripada kasus pendinginan karena osilasi menuntut langkah yang rapat. Kurva posisi berosilasi teredam dari 2 menuju nol dengan periode sekitar 3.1 detik yang sesuai $\sqrt{4 - 0.04}$ , sedangkan kurva kecepatan memimpin seperempat periode di depannya. Pemisahan y(:, 1) dan y(:, 2) pada plot menjadi kebiasaan yang perlu dipegang, karena plot(t, y) sekalipun menggambar kedua kolom tanpa label yang menjelaskan mana posisi dan mana kecepatan.

Mengatur toleransi dengan odeset#

Pemecah menghentikan pembaruan langkah ketika error lokal diperkirakan masuk dalam toleransi, dan kedua ambangnya diatur oleh odeset. Opsi RelTol mengatur error relatif terhadap besar solusi, AbsTol mengatur batas mutlak yang berlaku ketika solusi mendekati nol, dan nilai bakunya berturut-turut 1e-3 dan 1e-6. Ketentuan yang sama berlaku pada kedua implementasi, sehingga memperketat toleransi selalu dibayar dengan langkah yang lebih banyak, dan jumlah langkah adalah ukuran biaya yang dapat dicetak secara deterministik.

matlab Dua toleransi menghasilkan jumlah langkah yang berbeda

k = 0.08;
T_aset = 20;
pendinginan = @(t, T) -k*(T - T_aset);
opsi_longgar = odeset('RelTol', 1e-3, 'AbsTol', 1e-6);
opsi_ketat = odeset('RelTol', 1e-9, 'AbsTol', 1e-12);
[t1, T1] = ode45(pendinginan, [0 60], 90, opsi_longgar);
[t2, T2] = ode45(pendinginan, [0 60], 90, opsi_ketat);
fprintf('RelTol 1e-3 : %d langkah\n', numel(t1));
fprintf('RelTol 1e-9 : %d langkah\n', numel(t2));

text Keluaran

RelTol 1e-3 : 81 langkah
RelTol 1e-9 : 309 langkah

Memperketat toleransi dari $10^{-3}$ ke $10^{-9}$ melipatgandakan jumlah langkah hampir empat kali, dari 81 menjadi 309. Kedua solusi menggambar kurva yang tidak dapat dibedakan mata pada kasus sederhana ini. Pemilihan toleransi karenanya menjadi kompromi antara biaya dan ketelitian yang baru terasa pada simulasi panjang atau sistem yang sensitif. Kebiasaan yang sehat adalah memulai dari nilai baku, memeriksa apakah hasil berubah ketika toleransi diperketat satu orde, lalu berhenti pada toleransi yang hasilnya sudah stabil.

Kesalahan umum#

Menukar urutan keluaran menjadi [y, t]. Keduanya sama-sama matriks bernomor sehingga penugasan salah urut tidak memicu error, tetapi kolom pertama y kemudian dibaca sebagai waktu. Kebiasaan pencegahnya adalah membaca pasangan sebagai waktu dahulu baru solusi, sama seperti urutan pada rumus $y(t)$ .
Mengubah variabel yang ditangkap fungsi anonim. Nilai k dan T_aset disalin ketika handle pendinginan dibuat, sehingga menugaskan ulang k di workspace tidak mengubah model yang sudah berjalan. Buat ulang handle setelah parameter berubah, atau tulis parameter sebagai argumen model.
Menulis interval waktu terbalik tanpa sengaja. Rentang [t1 t0] yang menurun tetap sah dan memicu penyelesaian mundur dari t1, sehingga hasilnya berupa solusi pada waktu lampau alih-alih error. Periksa urutan batas ketika solusi yang diharapkan tidak muncul pada rentang waktunya.

Ringkasan#

Masalah nilai awal $y' = f(t, y)$ dengan $y(t_0) = y_0$ diselesaikan ode45 yang menerima model sebagai fungsi atau handle, rentang waktu, dan nilai awal kolom.
Keluaran [t, y] menempatkan waktu pada vektor kolom dan satu variabel keadaan pada satu kolom y, sehingga jumlah barisnya menyatakan jumlah langkah.
Persamaan orde dua diturunkan menjadi sistem dua orde satu dengan kecepatan sebagai keadaan kedua, dan model sistem ditulis paling rapi sebagai berkas fungsi tersendiri.
odeset mengatur RelTol dan AbsTol, dan memperketat toleransi meningkatkan jumlah langkah yang tercetak sebagai ukuran biaya deterministik.

Referensi#

ode45 — MathWorks — Rujukan resmi pemecah Runge-Kutta orde empat-lima beserta bentuk pemanggilannya.
odeset — MathWorks — Rujukan struktur opsi pemecah termasuk toleransi relatif dan absolut.