其他分享
首页 > 其他分享> > 数模-微分方程(人口预测之马尔萨斯模型和阻滞增长模型)

数模-微分方程(人口预测之马尔萨斯模型和阻滞增长模型)

作者:互联网

模型:

image
image

代码:

%% 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,'-')

image

image

例题

image

求解过程

image
image
image
image

代码

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)  % 预测的数值

image

标签:人口预测,xm,模型,拟合,数模,exp,year,x0,population
来源: https://www.cnblogs.com/jgg54335/p/15181619.html