GSL의 ODE 해결사를 사용하여 "MaxSteps"에 해당하는 것은 무엇입니까?

Question

Mathematica을 사용하여 생성 된 ODE 해석기를 GSL으로 재현하고 싶습니다. nbins는 해석 및 _h받는 주어진 방정식의 개수

int run_simulation() { 
    gsl_odeiv_evolve* e = gsl_odeiv_evolve_alloc(nbins); 
    gsl_odeiv_control* c = gsl_odeiv_control_y_new(1e-17, 0); 
    gsl_odeiv_step* s = gsl_odeiv_step_alloc(gsl_odeiv_step_rkf45, nbins); 
    gsl_odeiv_system sys = {function, NULL, nbins, this }; 
    while (_t < _tmax) { //convergence check here 
     int status = gsl_odeiv_evolve_apply(e, c, s, &sys, &_t, _tmax, &_h, y); 
     if (status != GSL_SUCCESS) { return status; } 
    } 
    return 0; 
} 

: GSL하여 정확한 상당려고

result[r_] := NDSolve[{ 
    s'[t] == theta - (mu*s[t]) - ((betaA1*IA1[t] + betaA2*IA2[t] + betaB1*IB1[t] + betaB2*IB2[t]) + 
            (betaA1T*TA1[t] + betaA2T*TA2[t] + betaB1T*TB1[t] + betaB2T*TB2[t])) * s[t] - 
           ((gammaA1*IA1[t] + gammaA2*IA2[t] + gammaB1*IB1[t] + gammaB2*IB2[t]) + 
            (gammaA1T*TA1[t] + gammaA2T*TA2[t] + gammaB1T*TB1[t] + gammaB2T*TB2[t])), 

//... Some other equations 



s[0] = sinit,IA1[0] = IA1init,IA2[0] = IA2init, 
IB1[0] = IB1init,IB2[0] = IB2init,TA1[0] = TA1init, 
TA2[0] = TA2init,TB1[0] = TB1init,TB2[0] = TB2init}, 
{s,IA1,IA2,IB1,IB2,TA1,TA2,TB1,TB2},{t,0,tmax}, 
MaxSteps->100000, StartingStepSize->0.1, Method->{"ExplicitRungeKutta"}]; 

: 여기

는 NDSolve 사용 티카 코드 현재 스텝 크기.

방정식 자체는 제공하지 않지만, 계단 숫자를 제한하는 방법 (Mathematica에서 MaxSteps->100000과 같이)을 찾은 유일한 방법은 gsl_odeiv_control_y_new 제어 기능의 첫 번째 인수를 수정하는 것입니다. 여기서 1e-17은 나에게 140000 개의 스텝을 제공합니다 ...

GSL의 ODE 해결사가 주어진 최대 단계 수를 사용하도록 강제하는 방법을 아는 사람이 있습니까? 당신이 아마 이해했듯이, 내가 그 두 도구들 사이에서 실제로 비교할 수있는 결과를 갖는 것이 중요합니다.

도움 주셔서 감사합니다.

Answer 1

MaxSteps은 RK (Runge Kutta)가 멈추었을 때만 중요하므로 결과적으로 시스템이 적절하게 발전하지 못합니다. 원하는 단계 수나 필요한 정확도가 수정되지는 않습니다. 물론 정확도가 높으면 단계 크기가 작아지기 때문에 일정한 간격으로 더 많은 단계가 필요합니다. 그러나 요점은 RK가 멈추거나 이상하지 않은 별난 시스템이 없거나 (이 경우 Mathematica 오류 메시지가 보이지 않는 경우) 또는 maxsteps를 작게 설정하면 MaxSteps가 mathematica를 적절하게 비교하는 데 도움이되지 않는다는 것입니다 및 GSL.

적절한 비교를하려면 두 프로그램에서 동일한 정확도 요구 사항 및 제어 기능을 설정해야합니다. 실제로 표준 옵션 외에도 API gsl_odeiv2_control_alloc 및 gsl_odeiv2_control_hadjust 기능을 사용하여 GSL의 임의 제어 기능을 설정할 수 있습니다. 매쓰 매 티카 코드에서 정확히 멈추는 조건이 무엇인지 확인해야합니다.

또 다른 옵션은 두 프로그램 모두에서 비 적응 고정 스텝 RK를 사용하는 것입니다 (gsl에서는 gsl_odeiv2_driver_apply_fixed_step을 호출하여 수정 단계로 시스템을 호출 할 수 있습니다).

마지막으로. 1e-17은 미친 상대 정확도 요구 인 것 같습니다. 반올림 오류는 일반적으로 RK가이 정도의 정확도에 도달하는 것을 허용하지 않는다는 것을 기억하십시오. 실제 반올림 오류는 RK가 멈추거나 Mathematica/GSL이 서로 동의하지 못하게 만드는 것 중 하나입니다. 정확도를 1e-10보다 높게 설정해야합니다.