2014-07-18 2 views
1

나는 이것을 research paper에 따르려고합니다. 20 페이지의 그림 7에있는 솔루션 그래프를 복제하려고합니다. 그림 7의 스크린 샷이 있습니다. enter image description hereMATLAB에서 비선형 미분 방정식의 시스템 해법을 그리는 방법은 무엇입니까?

먼저 왼쪽 그림을 다시 만들고 싶습니다. 문제의 시스템은 내가 가지고있는 것입니다 dX. 여기가 m-파일에있는 것입니다 :

function dX = CompetitionModel(t,X)  
    bs = 8*10^(-3); 
    bl = 4*10^(-3); 
    bh = 6.4*10^(-3); 
    N = bs + bl + bh; 
    K = 10^8; 
    m1 = 2*10^(-5); 
    m2 = 9*10^(-9); 
    p = 5*10^(-13); 
    I = 10^(-3); 
    T = 0; 
    a = 0; 
    dX = [X(1) * (bs * (1 - N/K) - I - T - m1) - p * X(1) * (X(2) + X(3)); 
     X(2) * (bl * (1 - N/K) - I - a*T - m2) + m1 * X(1) + p * X(2) * (X(1) - X(3)); 
     X(3) * (bh * (1 - N/K) - I - a*T) + m2 * X(2) + p * X(3) * (X(1) + X(2))]; 
end 

ode45 구문이 있습니다 [T,Y] = solver(odefun,tspan,y0). 나는 내가 게시 한 그림에서 tspan을 얻는다. 내 초기 조건은 : S0 = 10^4; Rl0 = 0; Rh0 = 0, 그래서 이것은 내가 가지고있는 것입니다 y0. 나는 명령 창에 다음을 입력 :

>>[t,X1] = ode45('CompetitionModel', [0,45000], [10^4, 0, 0]); 
>>[t,X2] = ode45('CompetitionModel', [0,45000], [10^4, 0, 0]); 
>>[t,X3] = ode45('CompetitionModel', [0,45000], [10^4, 0, 0]); 

MATLAB 지난 30 분 동안 바빴다 내 노트북이 정말 섹시지기 시작한다. 따라서 작업이 끝날 때까지 플로팅을 할 수 없으며 코드에 오류가 있는지 알 수 없습니다. 시스템의 솔루션을 얻을 수있는 더 좋은 방법이 있는지 궁금합니다. dX.

답변

3

나는 종이에 대한 ODE를 확인했는데 하나의 실수를 발견했다. 논문 9면 마지막 단락에 따르면 N = S + Rl + Rh. 수정 된 코드 :

function dX = CompetitionModel(~,X)  
    bs = 8e-3; 
    bl = 4e-3; 
    bh = 6.4e-3; 
    N = sum(X); % this line was incorrect 
    K = 1e8; 
    m1 = 2e-5; 
    m2 = 9e-9; 
    p = 5e-13; 
    I = 1e-3; 
    T = 0; 
    a = 0; 
    dX = [X(1) * (bs * (1 - N/K) - I - T - m1) - p * X(1) * (X(2) + X(3)); 
     X(2) * (bl * (1 - N/K) - I - a*T - m2) + m1 * X(1) + p * X(2) * (X(1) - X(3)); 
     X(3) * (bh * (1 - N/K) - I - a*T) + m2 * X(2) + p * X(3) * (X(1) + X(2))]; 

두 가지 구문상의 변화에 ​​유의하십시오. 먼저, ODE에서 시간 값을 실제로 사용하지 않기 때문에 함수 정의에서 사용하지 않은 입력에 대한 보조자로 사용한 t 대신 ~을 남겨 둘 수 있습니다. 둘째로 8*10^(-3) 대신에 8e-3이라는 표기법을 사용할 수 있습니다. 그들은 같은 것으로 평가하지만, 전자는 조금 더 깨끗해 보입니다.

무 처리의 플롯, 즉 T = 0을하기에 나타낸다.

No Treatment Population

내가 먼저 MATLAB에서 뻣뻣한 ODE 솔버들과 함께 당신의 방정식을 해결하기 위해 노력하여 ODE의 문제를 알아 냈어. 자세한 내용은 MATLAB ODE solver documentation을 참조하십시오. 기본적으로 나는 ode15sode23s으로 해결했고 솔루션은 불안정했다 (인구는 무한대로 떨어졌다). 이러한 다른 솔버는 솔루션이 걸려있을 때 명심해야 할 좋은 도구입니다. 때로는 다른 사람들 중 한 사람이 일하게 될 것이고, 당신이 원하는 것을 당신에게 주거나 다른 곳에 다른 문제가 있음을 나타낼 것입니다.

참고 : "위상 인물"을 올바르게 사용하지 않았 음을 확신합니다. 주어진 초기 조건 세트로 시간에 관한 ODE의 해를 찾고있을뿐입니다. phase portrait은 시스템의 상태 (여기서 세 개의 서로 다른 인구수)가 초기 조건의 여러 세트를 고려하여 어떻게 진화하는지 살펴 봅니다. 이 솔루션의 시간 의존성에 대해서는 전혀 알려주지 않았고, 다른 솔루션과 비교하여 어떻게 진화했는지는 알 수 없습니다.

+0

실수를 시정 해 주셔서 감사합니다. DE의 시스템 솔루션의 플롯으로 "위상 초상화"를 이해합니다. 너는 내 혼란을 해결했다. 시간 내 주셔서 감사합니다. –

관련 문제