Logistic模型的分岔图（MATLAB实现）

clc, clearx = 0.5i  =  1for mu=0:0.00001:2    x = 1-mu*x*x     派信   迅羡做y(i) = x     亩衡   i = i+1endplot(0:0.00001:2,y,'.')hold onxlabel('\mu')ylabel('x')虫口模型分岔图

怎么利用人口数据和非线性最小二乘法估计y=Xm/[1+(Xm/X0-1)*exp(-r*t)]式中r,x0,xm的值？其求解液陪山方法如下：

1、将【1790年乱御～2000年】的年份按序列来表示，闹中并赋值给t。即

t=1:22%1——代表1790年，2——代表1800年，。。。，22——代表2000年

2、将人口数赋值给y

3、自定义logistic模型函数，即

fun=inline('a(1)./(1+(a(1)/a(2)-1)*exp(-a(3).*t))','a','t')

4、使用nlinfit非线性最小二乘法函数，拟合系数Xm、X0、r，即

a= nlinfit(t,y,fun,b)

5、使用plot绘图函数，绘制拟合前后对比图

6、根据上述方法，编写程序并执行，可以得到如下结果。

【亩数1】初始值是这么定的：

>>[c,d]=solve('c/(1+(c/60.2-1)*exp(-5*d))=67.2','c/(1+(c/60.2-1)*exp(-20*d))=90.9','c','d')

c =

62.457687920920008140916437719297-1.6291275002135204129081793780683*i

62.457687920920008140916437719297+1.6291275002135204129081793780683*i

241.95986840768380162085441918688

d =

-.95675170588036288272009149214033e-1+.43712460153973020393407947582405*i

-.95675170588036288272009149214033e-1-.43712460153973020393407947582405*i

.29854891132649696390901254874996e-1

【2】fun的函数形式就是logistic函数。

【3】b1 =

180.9871 0.0336

可以改成这样，效果一样：

%==============

[c,d]=solve('c/(1+(c/60.2-1)*exp(-5*d))=67.2','c/(1+(c/60.2-1)*exp(-20*d))=90.9','c','d') %求初始参数

t=0:51%令1954年为初始年

x=[60.2 61.5 62.8 64.6 66 67.2 66.2 65.9 67.3 69.1 70.4 72.5 74.5 76.3 78.5 80.7 83 85.2 87.1 89.2 90.9 92.4 93.7 95 96.259 97.5 98.705 100.1 101.654 103.008 104.357 105.851 107.5 109.3 111.026 112.704 114.333 115.823 117.171 118.517 119.85 121.121 122.389 123.626 124.761 125.786 126.743 127.627 128.453 129.227 129.988 130.756]

b0=[ eval(c(3)), eval(d(3))]%初始参数值

fun=inline('b(1)./(1+(b(1)/60.2-1)*exp(-b(2)*t))','b','t')

[b1,r1,j1]=nlinfit(t,x,fun,b0)

y1=fun(b1,t)%非线性拟合的方程

plot(t+1954,x,'*',t+1954,y1) %对原始数据与曲线拟合后的值作图

grid

R1=r1.^2

R2=(x-mean(x)).^2

R=1-R1/R2 %可决虚信系数差耐轮

W=sum(abs(r1)) %残差绝对值之

%==============

---