【数模】【matlab】微分方程

人口增长模型

指数增长模型（Malthus人口模型）

模型假设

符号说明

模型建立

 t时刻人口总量为x(t)，对一个地区或国家来说，x(t)数量很大，视为一个可微函数。

参数估计

--------------------线性最小二乘

这类模型的目的是为了预报预测，预测准不准确不是模型本身可以完全决定的，需要用实际数据检验。所以，在数据拟合时，一般采用前2/3的数据来拟合参数（自找导师监督学习），用后1/3的数据来检验模型的可用性。如果数据本身较少，也说明可获取的信息太少，模型的可信度本来就低。 

 误差分析

优缺点

代码

clear
A=xlsread('D:\数据\美国人口.xls');
A=A';
t=A(:,1);x=A(:,2);    //注意行和列的赋值方式区别
t1=t(1:12);y=log(x(1:12));     //t1和y准备做线性回归
[b,bt,r,rt,st]=regress(y,[ones(length(y),1),t1]);     //
x0=exp(b(1));        //b（1）求出来是lnx0,就是第一个常数的大小
r=b(2);     //第二个常数
y1=x0*exp(r*t1);    //变换到原来的形式
plot(t1,x(1:12),'*',t1,y1,'-')

• 函数讲解   regress   多元线性回归

`regress` 是 MATLAB 中用于执行线性回归分析的函数。它的语法如下：

b------ returns a vector b of coefficient estimates for a multiple linear regression of the responses in vector y on the predictors in matrix x. The matrix x must include a column of ones. b是一个列向量返回的是在x和y这组数据下得到的多元线性回归的系数估计向量b，矩阵x必须包含一列1。至于x为什么需要包含一列1见下面对x的分析。

bint------returns a matrix bint of 95% confidence intervals for the coefficient estimates. 返回系数估计值的95%置信区间的矩阵绑定。bint是一个多行两列的矩阵，返回系数估计值的95%置信区间。对于置信区间的理解简单点说就是可信度，学过概率论就会明白何为置信区间。bint的每一行有两个值，可看成一个区间的上下界，它的每一行这个区间就是参数b中每一行的那个参数的95%的置信区间。

r------returns an additional vector r of residuals. r是residual的简写，就是残差的意思。r返回一个列向量，r中的每一个值就是真实的数据y减去预测的y值得到的，称之为残差。

rint------returns a matrix rint of intervals that can be used to diagnose outliers. 返回的也是一个多行两列的矩阵，它的理解可同bint，只不过此时判断的是r的可信区间。

stats------returns a vector stats that contains the R2   越接近1，回归方程越显著
 

y------是一个列向量，是已知数据。

x------是一个矩阵，它的第一列全为1，有没有全为1的这一列决定回归方程是否含有拟合的常数参数，也就是y=ax21+bx1+cx2+d中的d这个参数是否有被拟合出来。在程序中对应于全1的那一列一般表示为:ones(size(y)) ,其中x1和x2是多元线性回归方程中的自变量，是已知数据

• ones(length(y), 1) 的作用是创建一个全为 1 的列向量，用作解释变量矩阵的第一列。它在线性回归分析中的作用是添加一个常数项，通常表示截距。

通常，在线性回归中，我们将解释变量矩阵的第一列设置为全为 1 的列向量，以便模型考虑到截距项。这样可以确保即使当解释变量为0时，模型也能有非零的预测值。

regress详细请看

http://t.csdnimg.cn/WebVv

阻滞增长人口模型（Logistic模型）

 模型假设

 

模型略

主要就是增长率可改变，人口增长呈s型曲线

代码

function x=Logistcfunc(a,t)
x=a(1)./(1+(a(1)/a(2)-1)*exp(-a(3)*(t-1790)/10));



clear
A = xlsread('D:\数据\美国人口.xls');
A = A';
t = A(:, 1);
x = A(:, 2);
t1 = t(1:21);
x1 = x(1:21);
[a, re] = lsqcurvefit(@Logistcfunc, [300, 3, 0.02], t1, x1);
xt = a(1)./(1+(a(1)/a(2)-1)*exp(-a(3)*(t-1790)/10));
plot(t, xt, '-', t, x, '*')

 函数讲解  lsqcurvefit

格式: 

x = lsqcurvefit(fun,x0,xdata,ydata)

f: 符号函数句柄,如果是以m文件的形式调用的时候,别忘记加@.这里需要注意,f函数的返回值是和y匹对的,即拟合参数的标准是(f-y)^2取最小值,具体看下面的例子

x0:最开始预估的值(预拟合的未知参数的估计值)。如上面的问题如果我们预估A为1,B为2,C为3,则x0=[1 2 3]

xdata:我们已经获知的x的值

ydata:我们已经获知的x对应的y的值

扩展形式

[x,resnorm,residual,exitflag,output] = lsqcurvefit(fun,x0,xdata,ydata,lb,ub,options);

• x: 最优解得到的参数值。
• resnorm: 残差的平方和，即残差的 2-范数的平方。
• residual: 在最优解 x 处的残差，即拟合函数 fun(x,xdata) 减去 ydata 的值。
• exitflag: 描述算法终止条件的标志。
• output: 包含优化过程信息的结构体。

http://t.csdnimg.cn/sHPdt

案例2 Volterra食饵-捕食者模型

函数脚本：

function xt=shier(t,x)
k1=1;k2=0.5;b=0.1;c=0.02;
xt=[x(1).*(k1-b*x(2));x(2).*(-k2+c*x(1))];

clear
ts=0:0.1:15;x0=[25,2];
[t,x]=ode45('shier',ts,x0);
subplot(1,2,1)   // 定位第一个图
plot(t,x),grid,gtext('x(t)'),gtext('y(t)');   //添加网格，gtext添加标签
subplot(1,2,2)
plot(x(:,1),x(:,2)),grid

函数讲解

ode45() 

[t, Xt] = ode45(odefun, tspan, X0)

[t,y] = ode45(odefun,tspan,y0)（其中 tspan = [t0 tf]）求微分方程组 y′=f(t,y) 从 t0 到 tf 的积分，初始条件为 y0。解数组 y 中的每一行都与列向量 t 中返回的值相对应。

一般只允许输入；两个参数

odefun 是函数句柄，可以是函数文件名，匿名函数句柄或内联函数名
tspan 是区间 [t0 tfinal] 或者一系列散点[t0,t1,…,tf] 
X0 是初始值向量  
t 返回列向量的时间点
Xt 返回对应T的求解列向量

参考文章：

【matlab中ode45函数使用的说明】_ZHAO__JW的博客-CSDN博客

【精选】matlab ode45的使用_ode45在matlab中的用法-CSDN博客

matlab 函数句柄详解_matlab函数句柄_晓林爱学习的博客-CSDN博客

【Matlab 控制】微分方程 ode45() 求解并绘制曲线_matlab ode45函数-CSDN博客