혼란 스 캐터링 시뮬레이션의 값이 기본 경우와 일치하지 않습니다.

Question

스택 오버플로에 대한 첫 번째 게시물은 부드럽게하십시오. I는 사용에 의해 해결된다 움직임

M(dv/dt)=-grad V(r),  dr/dt=v, 

방정식의 4 차원 시스템에 의해 설명되는 전위 V (R)에 질량 M의 입자의 X 위치, Y 평면을 수행하는 코드를 작성 Runge Kutta 4th order method, 여기서 r = (x, y) 및 v = (vx, vy); 이제 입자의 상태 X에 의해 정의되고, Y 및 E는 평면에서 에너지 인 벡터 V와 속도의 크기가

|v|=sqrt(2(E-V(r))/M) 

하여 주어진 양의 X 축이 이루는 각도 세타 상기 전위 V (R)이 지금 현재

V(r)=x^2y^2exp[-(x^2+y^2)], 

주어진다 I가 초기 값

x(0)=3, 
y(0)=0.3905, 
vx(0)=0, 
vy(0)=-sqrt(2*(E-V(x(0), y(0)))), 

만들어 코드이고 E = 0.260 * (1/EXP (2))

,617

// RK4 
    #include <iostream> 
    #include <cmath> 

    // constant global variables 
    const double M = 1.0; 
    const double DeltaT = 1.0; 

    // function declaration 
    double f0(double t, double y0, double y1, double y2, double y3); // derivative of y0 
    double f1(double t, double y0, double y1, double y2, double y3); // derivative of y1 
    double f2(double t, double y0, double y1, double y2, double y3); // derivative of y2 
    double f3(double t, double y0, double y1, double y2, double y3); // derivative of y3 
    void rk4(double t, double h, double &y0, double &y1, double &y2, double &y3); // method of runge kutta 4th order 
    double f(double y0, double y1); //function to use 

    int main(void) 
    { 
     double y0, y1, y2, y3, time, E, Em; 
     Em = (1.0/(exp(2.0))); 
     E = 0.260*Em; 
     y0 = 3.0; //x 
     y1 = 0.3905; //y 
     y2 = 0.0; //vx 
     y3 = -(std::sqrt((2.0*(E-f(3.0, 0.0)))/M)); //vy 
     for(time = 0.0; time <= 400.0; time += DeltaT) 
     { 
      std::cout << time << "\t\t" << y0 << "\t\t" << y1 << "\t\t" << y2 << "\t\t" << y3 << std::endl; 
      rk4(time, DeltaT, y0, y1, y2, y3); 
     } 


     return 0; 
    } 

    double f(double y0, double y1) 
    { 
     return y0*y0*y1*y1*(exp(-(y0*y0)-(y1*y1))); 
    } 

    double f0(double t, double y0, double y1, double y2, double y3) 
    { 
     return y2; 
    } 

    double f1(double t, double y0, double y1, double y2, double y3) 
    { 
     return y3; 
    } 


    double f2(double t, double y0, double y1, double y2, double y3) 
    { 
     return 2*y0*((y0*y0)-1)*(y1*y1)*(exp(-(y0*y0)-(y1*y1)))/M; 
    } 

    double f3(double t, double y0, double y1, double y2, double y3) 
    { 
     return 2*(y0*y0)*y1*((y1*y1)-1)*(exp(-(y0*y0)-(y1*y1)))/M; 
    } 


    void rk4(double t, double h, double &y0, double &y1, double &y2, double &y3) // method of runge kutta 4th order 
    { 
     double k10, k11, k12, k13, k20, k21, k22, k23, k30, k31, k32, k33, k40, k41, k42, k43; 
     k10 = h*f0(t, y0, y1, y2, y3); 
     k11 = h*f1(t, y0, y1, y2, y3); 
     k12 = h*f2(t, y0, y1, y2, y3); 
     k13 = h*f3(t, y0, y1, y2, y3); 
     k20 = h*f0(t+h/2, y0 + k10/2, y1 + k11/2, y2 + k12/2, y3 + k13/2); 
     k21 = h*f1(t+h/2, y0 + k10/2, y1 + k11/2, y2 + k12/2, y3 + k13/2); 
     k22 = h*f2(t+h/2, y0 + k10/2, y1 + k11/2, y2 + k12/2, y3 + k13/2); 
     k23 = h*f3(t+h/2, y0 + k10/2, y1 + k11/2, y2 + k12/2, y3 + k13/2); 
     k30 = h*f0(t+h/2, y0 + k20/2, y1 + k21/2, y2 + k22/2, y3 + k23/2); 
     k31 = h*f1(t+h/2, y0 + k20/2, y1 + k21/2, y2 + k22/2, y3 + k23/2); 
     k32 = h*f2(t+h/2, y0 + k20/2, y1 + k21/2, y2 + k22/2, y3 + k23/2); 
     k33 = h*f3(t+h/2, y0 + k20/2, y1 + k21/2, y2 + k22/2, y3 + k23/2); 
     k40 = h*f0(t + h, y0 + k30, y1 + k31, y2 + k32, y3 + k33); 
     k41 = h*f1(t + h, y0 + k30, y1 + k31, y2 + k32, y3 + k33); 
     k42 = h*f2(t + h, y0 + k30, y1 + k31, y2 + k32, y3 + k33); 
     k43 = h*f3(t + h, y0 + k30, y1 + k31, y2 + k32, y3 + k33); 


     y0 = y0 + (1.0/6.0)*(k10 + 2*k20 + 2*k30 + k40); 
     y1 = y1 + (1.0/6.0)*(k11 + 2*k21 + 2*k31 + k41); 
     y2 = y2 + (1.0/6.0)*(k12 + 2*k22 + 2*k32 + k42); 
     y3 = y3 + (1.0/6.0)*(k13 + 2*k23 + 2*k33 + k43); 

    } 

여기서 문제는 주어진 초기 조건으로 코드를 실행할 때, 값이이 문제

에 의해 주어진 경우에 따라에 예상되는 것과 일치하지 않는 것이 무엇인지 그래픽해야 지금

주어진 초기 조건과 같이, 생각 나는 방법의 올바른 구현을 가지고하지만 난이 코드를 실행하면 입자가 가능성에서 멀리 이동하기 때문에 그래프가 일치하지 않는 이유를 모르겠어 .

도움이 될 것입니다.

Answer 1

경로가 날카로운 방향으로 혼란스럽게 보입니다. 이를 위해서는 적응성 단계 크기가 필요하며 일부 단계 크기 제어를 구현해야합니다. 각 단계를 단계 길이의 절반 인 두 단계로 비교하거나 Fehlberg 또는 Dormand-Price와 같이 더 높은 차수의 방법이 내장 된 방법을 사용하여 비교하십시오.

---

더 즉각적인 오류 :

• 불필요한 매직 넘버
을 피하기 위해 Em V(1,1)로 정의

• 초기 위치는, 당신은 차트를 오른쪽으로 읽으면,

  y0 = 3.0; 
  y1 = -0.3905+k*0.0010; 
  

  와 k=-1,0,1, 빼기 기호에 유의하십시오.

• 초기 속도는 수평이며, 운동 에너지는 해당 위치의 잠재적 에너지를 보완하도록 계산됩니다.따라서 이러한 변화와 (python에서) 내가 얻을 적응 해결사 플롯으로

  y2 = v0 = -(std::sqrt((2.0*(E-V(y0, y1)))/M)); 
  y3 = v1 = 0.0; 
  

import numpy as np 
from scipy.integrate import odeint 
import matplotlib.pyplot as plt 

# capture the structure of the potential 
f = lambda  r :  r*np.exp(-r); 
df = lambda  r : (1-r)*np.exp(-r); 
V = lambda y1,y2 : f(y1*y1)*f(y2*y2); 

M= 1.0 
Em = V(1.0,1.0); 
E = 0.260*Em; 

def prime(t,y): 
    x1,x2,v1,v2 = y   
    dV_dx1 = 2*x1*df(x1*x1)*f(x2*x2); 
    dV_dx2 = 2*x2*df(x2*x2)*f(x1*x1); 
    return [ v1, v2, -dV_dx1/M, -dV_dx2/M ]; 

# prepare and draw the contour plot 
X1,X0=np.ogrid[-4:3:100j,-4:3:100j] 
plt.contour(X0.ravel(), X1.ravel(), V(X0,X1), Em*np.arange(0,1,0.1), colors='k', linewidths=0.3) 
# display grid and fix the coordinate ranges 
plt.grid();plt.autoscale(False) 

for k in range(-1,1+1): 

    x01 = 3.0; 
    x02 = b = -0.3905 + 0.0010*k; 
    v01 = -((E-V(x01,x02))*2.0/M)**0.5; 
    v02 = 0.0; 
    print "initial position (%.4f, %.4f), initial velocity (%.4f, %.4f)" % (x01, x02, v01, v02) 

    t0 = 0.0 
    tf = 50.0 
    tol = 1e-10 
    y0 = [ x01, x02, v01, v02 ] 
    t = np.linspace(t0,tf,501); y = odeint(lambda y,t: prime(t,y) , y0, t) 

    plt.plot(y[:,0], y[:,1], label="b=%.4f" % b, linewidth=2) 

plt.legend(loc='best') 
plt.show()