其他分享
首页 > 其他分享> > 数学建模姜启源5.1人口增长——学习记录

数学建模姜启源5.1人口增长——学习记录

作者:互联网

目录

指数增长模型(Malthus马尔萨斯人口模型)

方法一:线性最小二乘法 

 方法二:数值微分的中点公式

方法三:改进的指数增长模型

 logistic模型

方法一:线性最小二乘

方法二:非线性最小二乘拟合


指数增长模型(Malthus马尔萨斯人口模型)

方法一:线性最小二乘法 

clear
clc
t=1790:10:2000;
t=(t-1790)./10;
x=[3.9 5.3 7.2 9.6 12.9 17.1 23.2 31.4 38.6 50.2 62.9 76 92 105.7 122.8 131.7 150.7 179.3 203.2 226.5 248.7 281.4];
y=log(x);%y=ln(x)
coefficient=polyfit(t,y,1);  %polyfit线性最小二乘法拟合;第三个输入项n,表示用n次函数拟合曲线
%y=rt+a,y=lnx,a=lnx0
%coefficient输出拟合函数的参数p1,p2...pn(r,a)
r=coefficient(1);
x0=exp(coefficient(2));
%x(t)=x0e^rt
xt=x0*exp(r.*t);
plot(t,x,'o',t,xt,'-')%输出原数据图,拟合后函数图
error=sum((xt-x).^2);%计算误差

 误差=3.4742e+04

 方法二:数值微分的中点公式

clear
clc
t=1790:10:2000;
t=(t-1790)./10;
x=[3.9 5.3 7.2 9.6 12.9 17.1 23.2 31.4 38.6 50.2 62.9 76 92 105.7 122.8 131.7 150.7 179.3 203.2 226.5 248.7 281.4];
%数值微分的中点公式
n=length(x);
dx0=(-3*x(1)+4*x(2)-x(3))/(2*10);%左端点
dxn=(x(n-2)-4*x(n-1)+3*x(n))/(2*10);%右端点
dxk=zeros(1,n);
for i=1:n
    if i==1
       dxk(1,i)=dx0;
    elseif i==n
       dxk(1,i)=dxn; 
    else
        dxk(1,i)=(x(i+1)-x(i-1))/(2*10);
    end
end
rk=dxk./x;%计算rk
r=mean(rk);%r为平均增长率
x0=x(1);%第一年人口数
xt=x0*exp(r.*t*10);
plot(t,x,'o',t,xt,'-')%输出原数据图,拟合后函数图
xlabel('t')
ylabel('x')

误差=2.2048e+04

方法三:改进的指数增长模型,不再设定人口增长率r为不变常数,而是随着时间增加而减少。

clear
clc
t=1790:10:2000;
t=(t-1790)./10;
rt=[0.2949 0.3113 0.2986 0.2969 0.2907 0.3012 0.3082 0.2452 0.2435 0.242... 
    0.2051 0.1914 0.1614 0.1457 0.1059 0.1059 0.1579 0.1464 0.1161 0.1004 0.1104 0.1349];
x=[3.9 5.3 7.2 9.6 12.9 17.1 23.2 31.4 38.6 50.2 62.9 76 92 105.7 122.8 131.7 150.7 179.3 203.2 226.5 248.7 281.4];
%r(t)=r0-r1t
%%
coefficient=polyfit(t,rt,1);  %polyfit线性最小二乘法拟合;第三个输入项n,表示用n次函数拟合曲线
%coefficient——P(1)*X^N + P(2)*X^(N-1) +...+ P(N)*X + P(N+1)
r=polyval(coefficient,t);%拟合出的年增长率
r0=coefficient(2);
r1=-coefficient(1);
figure(1)
plot(t,rt,'o',t,r,'-')
xlabel('t')
ylabel('r')
title('美国人口10年增长率t-r散点图及拟合直线')
%%
x0=x(1);%第一年人口数
xt=x0*exp(r0.*t-r1*t.^2/2);
figure(2)
plot(t,x,'o',t,xt,'-')%输出原数据图,拟合后函数图
xlabel('t')
ylabel('x')
title('改进的指数增长模型')

误差=1.1311e+03 

 

 logistic模型

方法一:线性最小二乘

clear
clc
t=1790:10:2000;
t=(t-1790)./10;
x=[3.9 5.3 7.2 9.6 12.9 17.1 23.2 31.4 38.6 50.2 62.9 76 92 105.7 122.8 131.7 150.7 179.3 203.2 226.5 248.7 281.4];
%% 计算各参数
%数值微分的中点公式:求函数x(t)在各分点的导数近似值
n=length(x);
dxk=zeros(1,n);
dxk(1,1)=(-3*x(1)+4*x(2)-x(3))/(2*10);%左端点;
dxk(1,n)=(x(n-2)-4*x(n-1)+3*x(n))/(2*10);%右端点 
for i=2:n-1
    dxk(1,i)=(x(i+1)-x(i-1))/(2*10);%中间的点
end
y=dxk./x;
coefficient=polyfit(x,y,1);%线性最小二乘法拟合,(dx/dt)/x=-(r/xm)*x+r,系数由次数从大到小输出
r=coefficient(2);
xm=-coefficient(2)/coefficient(1);
x0=x(1);%第一年人口数
%% 画图
figure(1)
rt=[0.2949 0.3113 0.2986 0.2969 0.2907 0.3012 0.3082 0.2452 0.2435 0.242... 
    0.2051 0.1914 0.1614 0.1457 0.1059 0.1059 0.1579 0.1464 0.1161 0.1004 0.1104 0.1349];
xt=xm./(1+(xm/x0-1).*exp(-r*t*10));
rt1=polyval(coefficient,x)*10;
plot(x,rt,'o',x,rt1,'-')
xlabel('x')
title('美国人口增长率x-r散点图及拟合直线')

figure(2)
plot(t,x,'o',t,xt,'-')%输出原数据图,拟合后函数图
xlabel('t')
ylabel('x')
title('用方法一估计参数')
error=sum((xt-x).^2);%计算误差

方法二:非线性最小二乘拟合

clear
clc
t=1790:10:2000;
t=(t-1790)./10;
x=[3.9 5.3 7.2 9.6 12.9 17.1 23.2 31.4 38.6 50.2 62.9 76 92 105.7 122.8 131.7 150.7 179.3 203.2 226.5 248.7 281.4];
f=@(a,t)a(3)./(1+(a(3)/a(2)-1).*exp(-a(1)*t));
%r=a(1);x0=a(2);xm=a(3);
a0=[0.28 3.9 352];%采用比较好的初值,拟合会更精确,比如用方法一已经得到的值
[q,r] = lsqcurvefit(f,a0,t,x);%非线性最小二乘拟合
% 输入参数:
% fun:需要拟合的函数,假定有n个需要拟合的参数,那么 q = [q1,q2,...,qn]
% q_0:表示用户给定的一个起始点
% xdata:函数的自变量
% ydata:函数的因变量
% 输出参数:
% q:表示求解得到的最优参数
% r:表示最小二乘的目标函数值,即残差。
plot(t,x,'o',t,f(q,t),'-')
xlabel('t');ylabel('x')
title('用方法二估计参数')
disp(['误差=',num2str(r)])%输出误差


误差=458.2113

预测2010年,改进的指数增长模型误差:1.7%

                       logistic模型(方法二)误差:-3.8%

      

标签:coefficient,5.1,姜启源,10,1790,建模,dxk,拟合,xt
来源: https://blog.csdn.net/m0_59179752/article/details/122728714