模型：

程式碼：

%% Malthus模型（馬爾薩斯模型）
clear;clc
x = dsolve('Dx=r*x','x(0)=x0','t')    % x = dsolve('Dx=r*x','x(t0)=x0','t')
% x = x0*exp(r*t)
% 怎麼把上面這個式子中的x0和r替換成確定的值？
x0 = 100; 
r = 0.1;
subs(x)

% 初始人口為1000，畫出50年且增長率分別為0.5%，1%，1.5% 一直到5%的人口變化曲線
x0 = 1000;
for i = 1:10
    r = 0.005*i;
    xx = subs(x);
    fplot(xx,[0,50])   % fplot函式可以繪製表示式的圖形
    hold on
end


%% 阻滯增長模型（logistic模型）
clear;clc
x = dsolve('Dx=r*(1-x/xm)*x','x(t0)=x0','t')   % 化簡後和書上的結果一樣
% mupad  % 把計算出來的結果貼上過去可以得到直觀的表示式
% 高版本可以使用實時指令碼
t0 = 0;
x0 = 1000;
xm = 10000;
r = 0.05;
xx = subs(x);    %  10000/(exp(log(9) - t/20) + 1)
fplot(xx,[0,200])  

% 如果不會用上面的fplot函式怎麼辦？
t = 0:200;
x = 10000 ./ (exp(log(9) - t/20) + 1);
plot(t,x,'-')
 

例題

求解過程

程式碼

createFit.m

function [fitresult, gof] = createFit(year, population)
[xData, yData] = prepareCurveData( year, population );

% Set up fittype and options.
ft = fittype( 'xm/(1+(xm/3.9-1)*exp(-r*(t-1790)))', 'independent', 't', 'dependent', 'x' );
opts = fitoptions( 'Method', 'NonlinearLeastSquares' );
opts.Display = 'Off';
opts.StartPoint = [0.2 500];

% Fit model to data.
[fitresult, gof] = fit( xData, yData, ft, opts );

% Plot fit with data.
figure( 'Name', 'untitled fit 1' );
h = plot( fitresult, xData, yData );
legend( h, 'population vs. year', 'untitled fit 1', 'Location', 'NorthEast' );
% Label axes
xlabel year
ylabel population
grid on
 

code.m

clear;clc
year = 1790:10:2000;
population = [3.9,5.3,7.2,9.6,12.9,17.1,23.2,31.4,38.6,50.2,62.9,76.0,92.0,106.5,123.2,131.7,150.7,179.3,204.0,226.5,251.4,281.4];
cftool  % 擬合工具箱
% (1) X data 選擇 year
% (2) Y data 選擇 population
% (3) 擬合方式選擇：Custom Equation (自定義方程)
% (4) 修改下方的方框為：x = f(t) = xm/(1+(xm/3.9-1)*exp(-r*(t-1790)))
% (5) 左邊的result一欄最上面顯示：Fit computation did not converge:即沒有找到收斂解，右邊的擬合圖形也表明擬合結果不理想
% (6) 點選Fit Options，修改非線性最小二乘估計法擬合的初始值(StartPoint), r修改為0.02，xm修改為500
% (7) 此時左邊的result一覽得到了擬合結果：r = 0.02735, xm = 342.4
% (8) 依次點選擬合工具箱的選單欄最左邊的檔案—Generate Code(匯出程式碼到時候可以放在你的論文附錄)，可以得到一個未命名的指令碼檔案
% (9) 在這個開啟的指令碼中按快捷鍵Ctrl+S，將這個檔案儲存到當前資料夾。
% (10) 在現在這個檔案中呼叫這個函式得到引數的擬合值和預測的效果
[fitresult, gof] = createFit(year, population) 
t = 2001:2030;
xm = 342.4;   
r =  0.02735;
predictions = xm./(1+(xm./3.9-1).*exp(-r.*(t-1790)));  % 計算預測值（注意這裡要寫成點乘和點除）
figure(2)
plot(year,population,'o',t,predictions,'.')  % 繪製預測結果圖
disp(predictions)  % 預測的數值