1 차원 파를 시뮬레이션하려고합니다.

Question

저는 고조파 발진기 체인을 사용하여 1 차원 파의 간단한 시뮬레이션을하려고합니다. 교과서에 따른 - - 위치 (x[i]=0 평형 가정), i 번째 발진기 x[i] 대한 미분 방정식을 밝혀이 하나

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

은 (미분은을 WRT되고 시간) 그래서 나는 수치 적으로 다음과 같은 알고리즘으로 역학을 계산하려했습니다.

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); 
    } 

할 수 : 여기 osc[i] 애트리뷰트 loc (절대 위치) vel (속도) acc (가속도)와 bloc (평형 위치)와 오브젝트와 i 번째의 발진기 dt 시간 증분된다 알고리즘 here의 알고리즘을 참조하십시오. 클릭하면 6 번째 오실레이터에 충격이 가해집니다. 당신은 웨이브를 전혀 생성하지 않을뿐만 아니라 증가하는 총 에너지로 모션을 생성한다는 것을 볼 수 있습니다 (심지어 댐핑을 추가하더라도!). 내가 도대체 ​​뭘 잘못하고있는 겁니까?

Answer 1

모든 답변은 유용한 정보를 제공합니다. 당신이해야 할 것은 : 업데이트를 수행하기위한

• 사용 ja72's 관찰 - 당신이 (일관된 수와 동일한 반복 배치의 위치에 따라 노드에서 가속도를 계산하기 위해 필요한 즉, $를 함께 사용하지 마십시오 서로 다른 반복 단계에 속하는 금액이므로 같은 가속도 표현에서 x⁽k⁾ $ 및 $ x⁽k + 1⁾ $)
• tom10은 원래 위치를 잊어 버렸습니다. 상대 위치 만 고려하면됩니다. -이 웨이브 방정식의 동작은 다각형 선에 적용된 라플라시안 평활화 필터와 같습니다. - 이웃 한 두 개의 직접 연결된 이웃을 사용하여 점을 완화하고 이웃에 의해 결정된 세그먼트의 중간쪽으로 당겨받습니다.
• 마지막으로, 에너지를 보존하려는 경우 (예 : 물 표면이 진동하는 것을 막지 마십시오), 그것은 symplectic integrator를 사용할 가치가 있습니다. Symplectic/semi-implicit 오일러가합니다. 이러한 종류의 시뮬레이션을 위해서는 고차원 적분기가 필요하지 않지만 시간이 있다면 Verlet 또는 Ruth-Forest를 사용하는 것이 좋습니다. 정확도면에서 Symplectic 및 Order 2 또는 3을 사용하는 것이 좋습니다. Runge-Kutta는 암시 적 오일러처럼 시스템에서 에너지를 끌어 내 (느리게), 따라서 고유 한 감쇠를 얻게됩니다. Explicit Euler는 결국 당신을 위해 운명을 철자 할 것입니다 - 그것은 나쁜 선택입니다 - 항상 (뻣뻣하거나 손상되지 않은 시스템의 경우). 거기에 다른 통합의 톤이있다, 그래서 가서 실험.

P. 당신은 모든 입자에 동시에 영향을 주어야하기 때문에 진정한 암시 적 오일러를 구현하지 않았습니다. 이를 위해서는 투영 된 공액 그라데이션 방법을 사용하고 야 코비 행렬을 유도하고 입자 문자열에 대해 희소 선형 시스템을 풀어야합니다. 명시 적 또는 반 암시 적 방법을 사용하면 애니메이션에서 기대하는 "리더 따라 가기"동작을 얻을 수 있습니다.당신이 더 많은 것을 원하는 경우에

지금, 나는 구현 SciLab이 답 (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'); 

을

Answer 2

시간이 지남에 따라 늘어나는 진폭은 사용중인 간단한 오일러 통합의 결과 일 수 있습니다. 더 짧은 시간 간격을 사용하거나 더 나은 에너지 보존 특성을 갖는 symplectic 오일러 semi-implicit Euler으로 전환 해보십시오.

홀수 번식 행위 (파도가 매우 천천히 전파되는 곳)는 입자 질량에 비해 스프링 정수 (k)가 너무 낮을 수 있습니다.

또한 게시 한 방정식이 k * m이 아닌 k/m을 포함해야하므로 약간 잘못되었습니다. 즉, 방정식은

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

(c.f. this Wikipedia page)이어야합니다. 그러나 이는 전체 전파 속도에만 영향을줍니다.

Answer 3

당신은 두 가지 중요한 방법에 코드에서 잘못 초기 방정식을 표현했습니다

첫째, 방정식은 상대 왜곡을 표현하고 있습니다; 즉, 방정식에서 x[i]==x[i-1]==x[i+1]이면 x[i]의 0과 관계없이 x"[i]=0과 관계가 없습니다. 코드에서 각 입자에 대한 평형으로부터 절대 거리를 측정합니다. 즉 방정식은 끝 부분에있는 단일 스프링과 같이 경계에서 오실레이터의 행을 고정하지만 일부는 "평형"지점에 고정 된 전체 스프링 세트를 시뮬레이션합니다.

둘째, osc[i].acc+=+k*(osc[i-1].loc-osc[i-1].bloc);과 같은 용어는 의미가 없습니다. 여기에 절대 입자 위치 옆에만 osc[i].acc을 설정하고 상대 위치에 있지 않습니다.

이 두 번째 문제는 아마 당신이 에너지를 얻는 이유 일 것입니다.하지만 이는 부정확 한 시뮬레이션을 수행하는 부작용 일 뿐이며 반드시 시뮬레이션에서 오류가 발생했다는 것을 암시하지는 않습니다. 현재 간단한 오일러 시뮬레이션에서 변경해야한다는 증거는 없습니다 (Nathan이 제안했듯이). 먼저 방정식을 수정 한 다음 필요에 따라 시뮬레이션 방법을 사용하여 멋지게 만듭니다.

대신, osc.loc[i]-osc.loc[i-1] 같은 용어로, 즉 상대 변위에 대한 방정식을 작성합니다.

코멘트 : 나는 거의 내가 대답 한 질문에 다른 사람의 답변에 댓글을 남기고하지하지만 나단과 ja72 두 단지의 주요 문제가되지 않습니다 것들에 초점을 맞추고있다. 처음으로 시뮬레이션 방정식이 정확하고, 오일러보다 더 매끄러운 접근 방식에 대해 걱정할 필요가있을 때 등을 업데이트하는 순서 등이있을 수 있습니다. 간단한, 선형, 1 차 방정식, 특히 약간의 감쇠가있는 경우, 시간 단계가 충분히 작 으면 직접 오일러 방법이 잘 작동하므로 먼저 이와 같이 작동 시키십시오. 하나 귀하의 가속

Answer 4

osc[i].acc 여기에 업데이트 된 위치 osc[i-1].loc과하지 않은 업데이트 아직 위치 osc[i+1].loc에 따라 달라집니다.

먼저 모든 노드에 대해 모든 가속도를 계산 한 다음 위치 및 속도를 별도 루프로 업데이트해야합니다.

주어진 시간, 위치 및 속도의 가속도를 반환하는 함수를 만드는 것이 가장 좋습니다.

// 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] 
} 

고차원 적분기 체계를 사용해야 할 수도 있습니다. 2 차 암시 적 Runge-Kutta부터 시작하십시오.

Answer 5

여기에 설명 된 수퍼 간단하고 빠른 방법이 있습니다. 한 줄의 코드와 한 개의 루프만으로 매우 현실적입니다.

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

그러나 나는 파도가 반전없이 끝 지점에서 반사되도록 경계 조건을 설정할 수 없습니다입니다.

내가 할 수있는 최선책은 자유로운 끝에 vel = 0을 의미하는 마지막 하나를 설정하는 것이지만, 그 다음 웨이브는 간단히 끝나고 내가 원하는대로 되돌릴 수 없습니다.

아무도 알아낼 수 없습니까?