Cố gắng mô phỏng sóng 1 chiều

Question

Tôi đang cố gắng tạo mô phỏng đơn giản về sóng 1 chiều với một chuỗi các bộ dao động điều hòa. Các phương trình vi phân cho vị trí x[i] của dao động thứ i (giả định trạng thái cân bằng tại x[i]=0) hóa ra là - theo sách giáo khoa - một trong những điều này:

m*x[i]''=-k(2*x[i]-x[i-1]-x[i+1])

(đạo hàm được WRT các thời gian) vì vậy tôi đã cố gắng tính toán động lực bằng thuật toán sau đây. Dưới đây osc[i] là dao động thứ i như một đối tượng với các thuộc tính loc (tuyệt đối vị trí), vel (vận tốc), acc (tăng tốc) và bloc (cân bằng vị trí), dt là tăng thời gian:

for every i: 
    osc[i].vel+=dt*osc[i].acc; 
    osc[i].loc+=dt*osc[i].vel; 
    osc[i].vel*=0.99;  
    osc[i].acc=-k*2*(osc[i].loc-osc[i].bloc); 

    if(i!=0){ 
     osc[i].acc+=+k*(osc[i-1].loc-osc[i-1].bloc); 
    } 
    if(i!=N-1){ 
     osc[i].acc+=+k*(osc[i+1].loc-osc[i+1].bloc); 
    } 

Bạn có thể xem thuật toán trong hành động here, chỉ cần nhấp để đưa xung vào bộ dao động thứ 6. Bạn có thể thấy rằng không chỉ nó không tạo ra một làn sóng nào cả, nhưng nó cũng tạo ra một chuyển động với tổng năng lượng ngày càng tăng (thậm chí nếu tôi thêm sự giảm chấn!). Tôi đang làm gì sai?

Answer 1

Tất cả các câu trả lời nhất định đều cung cấp một thỏa thuận hữu ích (tuyệt vời). gì bạn phải làm là:

• sử dụng ja72's quan sát để thực hiện các bản cập nhật - bạn cần phải chặt chẽ và tính gia tốc tại một nút dựa trên vị trí so với cùng lặp hàng loạt (tức là không pha trộn $ x⁽k⁾ $ và $ x⁽k + 1⁾ $ trong cùng một biểu thức gia tốc vì đây là các số thuộc các bước lặp khác nhau)
• quên vị trí ban đầu của bạn như tom10 đã nói - bạn chỉ cần xem xét các vị trí tương đối - hành vi của phương trình sóng này giống như bộ lọc làm trơn laplacian được áp dụng cho một đường đa giác - một điểm thoải mái khi sử dụng hai hàng xóm được kết nối trực tiếp, được kéo về giữa phân khúc được xác định bởi những người hàng xóm
• cuối cùng, nếu bạn muốn năng lượng của mình được bảo tồn (tức là không có mặt nước ngừng rung), nó có giá trị sử dụng một tích hợp symplectic. Euler Symplectic/bán ngầm sẽ làm. Bạn không cần các trình tích hợp bậc cao hơn cho loại mô phỏng này, nhưng nếu bạn có thời gian, hãy xem xét sử dụng Verlet hoặc Ruth-Forest vì chúng là cả hai đối tượng và theo thứ tự 2 hoặc 3 về độ chính xác. Runge-Kutta sẽ, giống như Euler ngầm, rút ​​năng lượng (chậm hơn nhiều) từ hệ thống, do đó bạn sẽ có được giảm xóc vốn có. Rõ ràng Euler sẽ đánh vần doom cho bạn cuối cùng - đó là một lựa chọn tồi - luôn luôn (đặc biệt là cho các hệ thống cứng hoặc không bị nén). Có rất nhiều nhà tích hợp khác, nên hãy tiếp tục và thử nghiệm.

P.S. bạn đã không triển khai đúng Euler ngầm vì nó đòi hỏi đồng thời ảnh hưởng đến tất cả các hạt. Đối với điều đó, bạn phải sử dụng phương pháp gradient liên hợp dự đoán, lấy được Jacobians và giải các hệ thống tuyến tính thưa thớt cho chuỗi hạt của bạn. Phương thức rõ ràng hoặc bán ngầm sẽ giúp bạn có được hành vi "theo dõi nhà lãnh đạo" mà bạn mong đợi cho một hoạt ảnh.

Bây giờ, nếu bạn muốn biết thêm, tôi thực hiện và kiểm tra câu trả lời này trong SciLab (quá lười biếng để chương trình nó trong C++):

n=50; 
for i=1:n 
    x(i) = 1; 
end 

dt = 0.02; 

k = 0.05; 
x(20) = 1.1; 
xold = x; 
v = 0 * x; 
plot(x, 'r'); 
plot(x*0,'r'); 
for iteration=1:10000 
    for i = 1:n 
     if (i>1 & i < n) then 
      acc = k*(0.5*(xold(i-1) + xold(i+1)) - xold(i)); 
      v(i) = v(i) + dt * acc; 
      x(i) = xold(i) + v(i) *dt; 
     end 
    end 
    if (iteration/500 == round(iteration/500)) then 
     plot(x - iteration/10000,'g'); 

    end 
    xold = x; 
end 
plot(x,'b'); 

Sự phát triển sóng được nhìn thấy trong các đồ thị xếp chồng lên nhau dưới đây:

Answer 2

Biên độ ngày càng tăng theo thời gian có thể là một tạo phẩm của tích hợp Euler đơn giản mà bạn đang sử dụng. Bạn có thể thử sử dụng dấu thời gian ngắn hơn hoặc thử chuyển sang semi-implicit Euler, còn gọi là Euler đối xứng, có tính năng bảo tồn năng lượng tốt hơn.

Đối với hành vi lan truyền lẻ (nơi sóng dường như lan truyền rất chậm), có thể là do bạn có quá thấp một hằng số lò xo (k) so với khối lượng hạt.

Ngoài ra, phương trình bạn đã đăng có vẻ hơi sai vì nó phải liên quan đến k/m chứ không phải k * m. Tức là, phương trình phải là

x[i]''=-k/m(2*x[i]-x[i-1]-x[i+1]) 

(c.f. this Wikipedia page). Tuy nhiên điều này chỉ ảnh hưởng đến tốc độ truyền tổng thể.

Answer 3

Bạn đã bày tỏ phương trình ban đầu không đúng cách trong mã của bạn theo hai cách quan trọng:

Đầu tiên, lưu ý rằng phương trình chỉ diễn tả tương biến dạng; tức là, trong phương trình, nếu x[i]==x[i-1]==x[i+1] thì x"[i]=0 bất kể khoảng cách x[i] là từ 0. Trong mã của bạn, bạn đang đo khoảng cách tuyệt đối từ trạng thái cân bằng cho mỗi hạt. Nghĩa là, phương trình sửa chữa hàng dao động chỉ ở ranh giới, giống như một lò xo duy nhất được giữ ở đầu, nhưng bạn đang mô phỏng một tập hợp toàn bộ các lò xo nhỏ, mỗi lò được cố định tại một số điểm cân bằng.

Thứ hai, các cụm từ của bạn như osc[i].acc+=+k*(osc[i-1].loc-osc[i-1].bloc); không có ý nghĩa nhiều. Ở đây bạn đang đặt osc[i].acc chỉ dựa trên tuyệt đối vị trí của hạt bên cạnh nó, chứ không phải trên tương đối vị trí giữa hai. Vấn đề thứ hai này có lẽ là lý do tại sao bạn đạt được năng lượng, nhưng đó chỉ là một tác dụng phụ của việc mô phỏng không chính xác, không nhất thiết ngụ ý rằng mô phỏng của bạn đang tạo ra lỗi. Hiện tại không có bằng chứng nào cho thấy bạn cần thay đổi từ mô phỏng Euler đơn giản (như Nathan đã đề xuất). Đầu tiên có được phương trình chính xác, và sau đó nhận được ưa thích với phương pháp mô phỏng nếu bạn cần.

Thay vào đó, hãy viết phương trình cho tương đối chuyển vị, ví dụ: với các cụm từ như osc.loc[i]-osc.loc[i-1].

nhận xét: Tôi hiếm khi nhận xét về câu trả lời của người khác cho câu hỏi tôi đã trả lời, nhưng cả Nathan và ja72 đều tập trung vào những thứ không phải là vấn đề chính. Đầu tiên là để phương trình mô phỏng của bạn chính xác, rồi, có thể lo lắng về các cách tiếp cận fancier hơn Euler và thứ tự cập nhật các thuật ngữ của bạn, v.v ... nếu bạn cần. Đối với một phương trình đơn giản, tuyến tính, đầu tiên, đặc biệt là với một chút giảm xóc, phương pháp Euler trực tiếp hoạt động tốt nếu bước thời gian đủ nhỏ, do đó, làm cho nó hoạt động như thế này trước tiên.

Answer 4

Đối với một tăng tốc của bạn osc[i].acc tại đây phụ thuộc vào vị trí được cập nhật osc[i-1].loc và vị trí chưa được cập nhật osc[i+1].loc.

Bạn cần phải tính toán tất cả các gia tốc đầu tiên cho tất cả các nút, và sau đó trong một vòng lặp riêng biệt cập nhật vị trí và vận tốc.

Cách tốt nhất là tạo một hàm trả về thời gian tăng tốc, vị trí và vận tốc.

// calculate accelerations 
// each spring has length L 
// work out deflections of previous and next spring 
for(i=1..N) 
{ 
    prev_pos = if(i>1, pos[i-1], 0) 
    next_pos = if(i<N, pos[i+1], i*L) 
    prev_del = (pos[i]-prev_pos)-L 
    next_del = (next_pos-pos[i])-L 
    acc[i] = -(k/m)*(prev_del-next_del) 
} 

// calculate next step, with semi-implicit Euler 
// step size is h 
for(i=1..N) 
{ 
    vel[i] = vel[i] + h*acc[i] 
    pos[i] = pos[i] + h*vel[i] 
} 

Bạn có thể cần sử dụng lược đồ trình tích hợp thứ tự cao hơn. Bắt đầu bằng lệnh thứ hai tiềm ẩn Runge-Kutta.

Answer 5

Có một cách siêu đơn giản và nhanh hơn để thực hiện việc này ở đây. chỉ một dòng mã và một vòng lặp và rất thực tế.

http://users.softlab.ece.ntua.gr/~ttsiod/wavePhysics.html#

tuy nhiên tôi không thể thiết lập các điều kiện biên như vậy mà các sóng dội từ các điểm kết thúc mà không đảo ngược.

điều tốt nhất tôi có thể làm là đặt cái cuối cùng thành cái trước, có nghĩa là vel = 0 ở cuối tự do, nhưng sau đó sóng chỉ kết thúc và không bị lật ngược trở lại như tôi mong muốn.

Có ai tìm ra điều này không?