2014-03-24 31 views
0

我試圖用Scilab(我是總新手)對Lotka-Volterra模型進行參數估計。當我嘗試運行腳本時,Scilab警告說不連貫的減法。我想我的問題與this topic相同,但是那裏的解決方案使用了Matlab函數。與Scilab的Lotka Volterra模型的參數估計

這裏是我的腳本:

// 1. Create Lotka Volterra function 

function [dY]=LotkaVolterra(t,X,c,n,m,e) 
    IngestC = c * X(1) * X(2) 
    GrowthP = n * X(1) 
    MortC = m * X(2) 
    dY(1) = GrowthP - IngestC 
    dY(2) = IngestC * e - MortC 
endfunction 

// 2. Define the Nonlinear Least Squares functions 

function f = Differences (x) 
    // Returns the difference between the simulated differential 
    // equation and the experimental data. 
    c = x(1) ;n = x(2);m = x(3);e = x(4);y0 = y_exp(1,:);t0 = 0 
    y_calc=ode(y0',t0,t,list(LotkaVolterra,c,n,m,e)) 
    diffmat = y_calc' - y_exp 
    f = diffmat(:) 
endfunction 

function val = L_Squares (x) 
    // Computes the sum of squares of the differences. 
    f = Differences (x) 
    val = sum(f.^2) 
endfunction 

// Experimental data 
t = [0:19]'; 
H=[20,20,20,12,28,58,75,75,88,61,75,88,69,32,13,21,30,2,153,148]; 
L=[30,45,49,40,21,8,6,5,10,20,33,34,30,21,14,8,4,4,14,38]; 
y_exp=[H',L']; 


// compute the model cost function 
function [f, g, ind] = modelCost (x, ind) 
    f = L_Squares (x) 
    g = derivative (L_Squares , x) 
endfunction 

// use of optim function with loops to avoid local minimum 
tic 
i=0 
fitminx=zeros(4,100); 
fitminy=zeros(1,100); 
for c=[0:0.1:1] 
    for n=[0:0.1:1] 
     for m=[0:0.1:1] 
      for e=[0:0.1:1] 
       i=i+1 
       x0 = [c;n;m;e] 
       [ fopt , xopt , gopt ] = optim (modelCost , x0) 
       fitminx(:,i)=xopt; 
       fitminy(:,i)=fopt; 
      end 
     end 
    end 
end 
[a,b]=min(fitminy) 
fitminx(:,a) 
toc 

的錯誤信息是:

lsoda-- at t (=r1), mxstep (=i1) steps 
needed before reaching tout 
     where i1 is :  500             
     where r1 is : 0.4145715729197D+01          
Attention : Le résultat est peut être inexact. 

!--error 9 
Soustraction incohérente. 
at line  4 of function Differences called by : 
at line  2 of function L_Squares called by : 
at line  16 of function %R_ called by : 
at line  15 of function %deriv1_ called by : 
at line  58 of function derivative called by : 
at line  3 of function modelCost called by : 
       [ fopt , xopt , gopt ] = optim (modelCost , x0) 

感謝你給的興趣和時間去我的問題(對不起,我的英語)

+0

請複製確切的錯誤信息。即使它是法語,這比描述它還好。 –

回答

0

問題出在

diffmat = y_calc' - y_exp 

添加以下代碼:

disp("Y_calc dimensions:"); 
disp(size(y_calc')); 

disp("y_exp dimensions:"); 
disp(size(y_exp));  

發現:

Y_calc dimensions:  
91. 2. 

y_exp dimensions:  
20. 2. 

我一無所知預期的行爲和預期的矩陣大小,但至少是錯誤的根本原因。

+0

對不起,我犯了一個錯誤重寫代碼:t = [0:19]'。我改變了話題。不管怎麼說,還是要謝謝你 – Wendac

1

重複我的回答here

的問題是,不知何故求解達到一個地步,它不能解決每個t的頌歌,並在某一點停止。因此您的y_calc的尺寸小於y_exp

如果這不是你的問題,它只是通過更改Differences功能的6線diffmat解決了

diffmat = y_calc' - y_exp(1:size(y_calc',1),:)