Oscillating one-dimensional systems

Studi kasus dan Terjemahan

Ref. Linge S, Langtangen HP, Programming for Computations - A Gentle Introduction to Numerical Simulations with Python

Terjemahan

4.3.1 Penurunan Model yang Sederhana

alt text

Banyak sistem keteknikan (engineering) berkaitan dengan osilasi, dan persamaan diferensial merupakan kunci utama untuk memahami, memprediksi, dan mengontrol osilasi. Kita mulai dengan model paling sederhana yang berkaitan dengan dinamika penting dari sistem osilasi. suatu benda dengan massa m melekat/dikaitkan pada pegas dan bergerak sepanjang garis tanpa gesekan, lihat Gambar 4.15 di samping untuk sketsa (rolling wheels menunjukkan “tidak ada gesekan”). Ketika pegas diregangkan (atau dikompresi), gaya pegas menarik (atau mendorong) bodi (penampang m) kembali dan bekerja "melawan" gerakan. Lebih tepatnya, misalkan x (t) adalah posisi bodi pada sumbu x, dimana bodi bergerak. Pegas tidak direntangkan ketika x= 0, sehingga gaya adalah nol, dan x= 0 karenanya posisi keseimbangan bodi. Gaya pegas adalah -kx, dimana k adalah konstanta yang diukur. Kami berasumsi bahwa tidak ada gaya lain (mis., Tidak ada gesekan). Hukum Newton ke-2 F=ma kemudian memiliki F=-kx dan a=x ̈ ,

yang dapat ditulis ulang sebagai:

dengan memperkenalkan ω=√(k/m)   (yang sangat umum).

Persamaan (4.42) adalah persamaan diferensial orde kedua, dan oleh karena itu kita memerlukan dua kondisi awal, satu pada posisi x(0) dan satu pada kecepatan x’(0). Di sini kita memilih bodi untuk berhenti, tetapi menjauh dari posisi setimbang:

Solusi tepat untuk Pers. (4.42) dengan kondisi awal ini adalah x(t)=X0 cosωT. Ini dapat dengan mudah diverifikasi dengan mensubsitusikan ke Pers. (4.42) dan memeriksa kondisi awal. Solusinya mengatakan bahwa sistem massa pegas berosilasi bolak-balik seperti yang dijelaskan oleh kurva kosinus.

Persamaan diferensial (4.42) muncul dalam banyak konteks lainnya. Contoh klasik adalah pendulum sederhana yang berosilasi bolak-balik. Buku-buku fisika berasal, dari hukum gerak kedua Newton, itu diperoleh:

dimana m adalah massa bodi di ujung pendulum dengan panjang L, g adalah percepatan gravitasi, dan ϴ merupakan sudut yang dibuat pendulum dengan vertikal. Mempertimbangkan sudut kecil ϴ, sin ϴ ≈ ϴ, dan kita dapatkan Pers. (4.42) dengan x = ϴ, ω=√(g/L) , x(0)=Θ, dan x’(0)=0, jika Θ merupakan sudut awal dan pendulum diam di t=0.

4.3.2 Solusi Numerik

Kita telah melihat metode numerik untuk mengendalikan turunan orde kedua, dan beberapa pilihan lainnya merupakan tambahan, akan tetapi kita mengetahu cara menyelesaikan persamaan turunan orde pertama dan bahkan sistem-sistem pada persamaan orde pertama. Dengan hanya sedikit, tetapi cukup umum, cara yang dapat kita tuliskan pada persamaan 4.42 sebagai sebuah sistem orde pertama dari 2 persamaan turunan. Kita memperkenalkan u=x dan v=x^'=u' sebagai 2 fungsi baru yang tidak diketahui. Dua persamaan yang sesuai muncul dari definisi v=u' dan persamaan asal (4.42):

(memperlihatkan bahwa kita dapat menggunakan u"=v') untuk menghilangkan turunan orde kedua dari hokum kedua newton). Selanjutnya kita dapat menerapkan metode forward euler untuk persamaan 4.43 dan 4.44, seperti yang sudah dilakukan pada section 4.2.2:

Sehingga menghasilkan skema komputasi sebagai berikut,

4.3.3 Memprogram Metode Numerik; Kasus Khusus

Pertama-tama, bahkan sebelum mencoba menjalankan program, Anda harus menghitung dua langkah dalam putaran waktu dengan kalkulator sehingga Anda memiliki beberapa hasil antara untuk dibandingkan. Menggunakan X0 = 2. Dt = 0: 157079632679, dan ω = 2, kita mendapatkan u1 = 2, v = -1,25663706, u2 = 1,80260791, dan v2 = 2,51327412. Perhitungan semacam itu menunjukkan bahwa program itu tampaknya benar. (Kemudian, kita dapat menggunakan nilai-nilai tersebut untuk membangun tes unit dan fungsi tes yang sesuai.)

Langkah selanjutnya adalah mengurangi delta t parameter diskritisasi dan melihat apakah hasilnya menjadi lebih akurat. Gambar 4.17 menunjukkan solusi numerik dan tepat untuk kasus delta t = P / 40; P / 160; P / 2000. Hasilnya jelas menjadi lebih baik, dan resolusi terakhir memberikan grafik yang tidak dapat dibedakan secara visual. Namun demikian, resolusi terakhir melibatkan 6000 interval komputasi secara total, yang dianggap cukup banyak. Namun, ini bukan masalah pada laptop modern, karena perhitungan hanya membutuhkan sepersekian detik.

Meskipun 2000 interval per periode osilasi tampaknya cukup untuk solusi numerik yang akurat, grafik kanan bawah pada Gambar 4.17 menunjukkan bahwa jika kita meningkatkan waktu simulasi, di sini hingga 20 periode, ada sedikit pertumbuhan amplitudo, yang menjadi signifikan dari waktu ke waktu. . Kesimpulannya adalah bahwa metode Forward Euler memiliki masalah mendasar dengan amplitudo yang tumbuh, dan bahwa diperlukan delta yang sangat kecil untuk mencapai hasil yang memuaskan. Semakin lama simulasi, semakin kecil Delta t. Sudah pasti saatnya untuk mencari metode numerik yang lebih efektif!

4.3.4 Sebuah Penyelesaian dari Metode Numerik

Dalam skema Forward Euler,

kita dapat mengganti un pada persamaan terakhir dengan nilai unC1 yang baru dihitung dari persamaan pertama:

Sebelum membenarkan perbaikan ini secara matematis, mari kita coba pada contoh sebelumnya. Hasilnya muncul pada Gambar 4.18. Kita melihat bahwa amplitudo tidak tumbuh, tetapi fase tidak sepenuhnya benar. Setelah 40 periode (Gbr. 4.18 kanan) kita melihat signifikan perbedaan antara solusi numerik dan tepat. Penurunan t menurun kesalahan. Misalnya, dengan 2000 interval per periode, kami hanya melihat fase kecil kesalahan bahkan setelah 50.000 periode (!). Kita dapat menyimpulkan bahwa perbaikan tersebut menghasilkan metode numerik yang sangat baik! Mari kita tafsirkan skema yang disesuaikan secara matematis. Pertama kami memesan (4,49) - (4,50) sedemikian rupa sehingga perbedaan pendekatan terhadap derivatif menjadi transparan: (4,51)

 (4,52) Kami menafsirkan (4,51) sebagai persamaan diferensial sampel pada titik mesh tn, karena kami memiliki vn di sisi kanan. Sisi kiri kemudian perbedaan maju atau Meneruskan perkiraan Euler ke turunan u0 , lihat Gambar 4.2. Di samping itu, kami menginterpretasikan (4,52) sebagai persamaan diferensial sampel pada titik mesh tnC1, karena kami miliki di sisi kanan. Dalam hal ini, perbedaan aproksimasi pada sisi kiri adalah perbedaan ke belakang,

Gambar 4.19 mengilustrasikan perbedaan mundur. Kesalahan dalam perbedaan mundur sebanding dengan t, sama seperti untuk perbedaan maju (tetapi konstanta proporsionalitas dalam istilah kesalahan memiliki tanda yang berbeda). Diskretisasi yang dihasilkan metode untuk (4,52) sering disebut sebagai skema Backward Euler. Untuk meringkas, gunakan perbedaan maju untuk persamaan pertama dan mundur Perbedaan untuk hasil persamaan kedua dalam metode yang jauh lebih baik daripada hanya menggunakan maju perbedaan dalam kedua persamaan. Cara standar untuk mengekspresikan skema ini dalam fisika adalah dengan mengubah urutan persamaan,

dan terapkan perbedaan maju ke (4,53) dan perbedaan mundur ke (4,54):

Artinya, pertama kecepatan v diperbarui dan kemudian posisi u, menggunakan kecepatan yang paling baru dihitung. Tidak ada perbedaan antara (4,55) - (4,56) dan (4,49) - (4,50) sehubungan dengan akurasi, jadi urutan persamaan diferensial asli tidak apa-apa. Skema (4.55) - (4.56) berada di bawah nama Semi-implisit Euler4 atau Euler-Cromer. Implementasi (4.55) - (4.56) ditemukan dalam file osc_EC.py. Inti dari kode itu seperti

4.3.7 Metode Runge-Kutta Orde 4

Metode Runge-Kutta Orde 4 adalah metode yang sering digunakan secara luas untuk menyelesaikan ODEs, karena menghasilkan data dengan tingkat akurasi yang tinggi bahkan dalam time step yang tidak terlalu kecil.

Gambar 4.24

Algoritma; Pertama-tama kita nyatakan algoritma 4-stage

Rumus 4.59

Dimana Rumus 4.60 Rumus 4.61 Rumus 4.62

Aplikasi; Kita bisa menjalankan simulasi seperti pada Figs. 4.16, 4.18, dan 4.21, untuk 40 periode. 10 periode terakhir ditunjukan melalui Fig. 2.24. Hasil yang ditunjukan terlihat impresif sebagaimana penggunaan metode Euler-Cromer.

Implementasi; Tingkatan dalam metode Runge-Kutta orde-4 bisa dengan mudah diimplementasikan sebagai modifikasi dari osc_Heun.py code. Sebagai alternatif, salah satu dapat menggunakan osc_odespy.py code dengan menyediakan argumen odespy_methods-[odespy.RK4] untuk membandingkan fungsi.

Derivasi; Derivasi dari metode Runge-Kutta orde-4 dapat disajikan dengan cara pedagogis yang menyatukan banyak elemen fundamental dari teknik diskritisasi numerik dan bisa menggambarkan banyak aspek “numerical thinking ”ketika membangun perkiraan metode solusi.

Kita mulai dengan mengintegrasikan general ODE (gambar rumus1) dari waktu ke waktu, mulai dari (gambar rumus2) sampai (gambar rumus 3),

(Persamaan1) Tujuan dari komputasi (persamaan2), ketika (persamaan3) pada saat ini lebih dikenal dengan nilai (u tanda miring). Tantangan mengintegralkan muncul ketika integrand mengandung u yang tidak diketahuai antara tn sampai tn+1.

Integral tersebut dapat diperkirakan dengan menggunakan Simpson’s rule yang telah terkenal (Persamaan 4) Permasalahan dengan persamaan ini adalah kita tidak mengetahui nilai dari (persamaan5) dan (pers6) karena hanya u^n yang tersedia dan hanya f^n yang dapat dihitung. Untuk melanjutkana, idenya dalah menggunakan berbagai perkiraan untuk (persamaan7) dan (pers8) berdasarkan penggunaan skema yang telah diketahui untuk ODE dalam interval (pers9) dan (pers10). Mari kita bagi persamaan integral menjadi empat suku. persamaan

4.3.9 ilustrasi redaman linier

Kami menganggap sistem rekayasa dengan pegas linier, s (u) = kx, dan peredam kental, di mana gaya peredaman adalah porpotional terhadap u ', f (u') = bu ', untuk beberapa konstanta b> 0. Pilihan ini dapat memodelkan sistem pegas vertikal di dalam mobil (tetapi insinyur sering suka menggambarkan sistem tersebut dengan massa bergerak horizontal seperti yang digambarkan pada Gambar 4.25). kita dapat memilih nilai-nilai sederhana untuk konstanta untuk mengilustrasikan efek dasar redaman (dan kegembiraan selanjutnya). Memilih osilasi sebagai fungsi u (t) = cos t sederhana dalam kasus undamped, kita dapat menetapkan m = 1, k = 1, b = 0,3, Uo = 1, Vo = 0. Fungsi berikut mengimplementasikan kasus ini:

Fungsi plot_u adalah kumpulan pernyataan plot untuk merencanakan u (t), atau bagian darinya. Gambar 4.27 menunjukkan efek dari istilah bu: kita memiliki osilasi dengan (perkiraan) periode 2π, seperti yang diharapkan, tetapi amplitudo teredam secara efisien.

Ceritakan tentang kebingungan dengan masalah berskala Daripada menetapkan b = 0,3 dan m - k = Uo = 1 sebagai nilai fisik yang cukup "tidak mungkin", akan lebih baik untuk skala persamaan mu + bu '+ ku = 0. Ini berarti bahwa kita memperkenalkan independen tanpa tergantung dan tidak tergantung variabel:

Di mana tc dan uc adalah ukuran karakteristik waktu dan perpindahan, masing-masing, sehingga ¯t dan ¯u memiliki ukuran tipikal mereka di sekitar kesatuan. Dalam masalah ini, kita dapat memilih uc = Uo dan tc = √ (m / k). Ini memberikan masalah skala (atau tanpa dimensi) berikut untuk kuantitas tak berdimensi ¯u (¯t):

Fakta yang mencolok adalah bahwa hanya ada satu parameter fisik dalam masalah ini: angka β yang tidak disebutkan. Memecahkan masalah ini terkait dengan penyelesaian masalah asli (dengan dimensi) dengan parameter m = k = Uo = 1 dan b = β. Namun, menyelesaikan masalah tak berdimensi lebih umum: jika kita memiliki solusiu (¯t; β), kita dapat menemukan solusi fisik dari berbagai masalah karena

Selama β diperbaiki, kita dapat menemukan u untuk Uo, k dan m dari rumus di atas! Dengan cara ini, simulasi yang memakan waktu dapat dilakukan hanya sekali, tetapi masih memberikan banyak solusi. Ini menunjukkan kekuatan bekerja dengan masalah skala atau tanpa dimensi.

4.3.10. Ilustrasi Redaman Linier Dengan Eksitasi Sinusoidal

Sekarang kita akan memperluas contoh sebelumnya untuk menambah beberapa gaya osilasi eksternal pada sistem: F (t) = Asin (wt). Mengendarai mobil di jalan dengan lonjakan sinusoidal mungkin memberikan eksitasi eksternal pada sistem pegas di mobil (w terkait dengan kecepatan mobil).

File:Coding4310.jpg

alt text

4.3.13 Metode finite diference; damping linier

Sebuah isu kunci adalah bagaimana untuk mengkonferensi skema dari daerah 4.3.12 ke persamaan diferensial dengan lebih banyak istilah. Kita mulai dengan kasus linear penempatan f (u') = bu', kemungkinan gaya per nonlinear s(u), dan sebuah gaya excitation F(t):

Kita harus cari perkiraan perbedaan yang tepat untuk u' di dalam bu'. Sebuah pilihan yang baik adalah perbedaan berpusat

Sampling persamaan pada titik tn,

dan memasukkan perkiraan perbedaan terhingga pada u" dan u' hasil dalam

dimana F" adalah notasi pendek untuk F(t). Persamaan (4.81) adalah linear dalam u^(n+1) tak diketahui kita dapat dengan mudah memecahkan untuk kuantitas ini: