数学建模之微分方程相关代码学习

Logistic模型求解

 1）非线性最小二乘估计

把第一个数据作为初始条件，利用余下的数据拟合参数Xm和r。

clc;
clear;
a=[1790:10:2000];%年份
%人口数据
b=[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,226.5,251.4,281.4];
b=b';%变成列向量
a=a';
a0=a(1);b0=b(1);%第一组数据作原始数据
fun=@(cs,td)cs(1)./(1+(cs(1)/b0-1)*exp(-cs(2)*(td-a0)));%cs(1)=Xm,cs(2)=r
cs=lsqcurvefit(fun,rand(2,1),a(2:end),b(2:end),zeros(2,1));%从第二组数据开始进行拟合
Y=fun(cs,[a;2010])%预测已知年代和2010的数据
c=[1790:10:2010];
plot(c,Y,'-r',a,b,'--b');
str=sprintf('Xm的拟合值为%f，r的拟合值为%f',cs(1),cs(2));
disp(str)

 

 2）线性最小二乘法

在这里十年一数据，故步长为10；

clc;
clear;
a=[1790:10:2000];%年份
%人口数据
b=[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,226.5,251.4,281.4];
b=b';%变成列向量
a=a';
a=[ones(21,1),-b(2:end)];
c=diff(b)./b(2:end)/10;%差分函数
cs=a\c;
r=cs(1);Xm=r/cs(2);
str=sprintf('Xm的拟合值为%f，r的拟合值为%f',Xm,r);
disp(str);

 再试试向前差分：

clc;
clear;
a=[1790:10:2000];%年份
%人口数据
b=[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,226.5,251.4,281.4];
b=b';%变成列向量
a=a';
a=[ones(21,1),-b(1:end-1)];
c=diff(b)./b(1:end-1)/10;%差分函数,往前差分
cs=a\c;
r=cs(1);Xm=r/cs(2);
str=sprintf('Xm的拟合值为%f，r的拟合值为%f',Xm,r);
disp(str);

 三种拟合方法，拟合同样的参数，由于方法不同，可能结果相差较大。

Matlab求微分方程的符号解

 微分方程的符号解

例一： （求解常微分方程的通解）

 

 

clc;
clear;
syms y(x)%定义符号变量
dsolve(x^2+y+(x-2*y)*diff(y)==0);
%diff()对于符号型变量为求导，对于数值型变量为差分
%一个“=”为赋值，==为等于；

例二：（求解常微分方程的初边值问题）

clc;
clear;
syms y(x)%定义符号变量
dy=diff(y);%一阶导
d2y=diff(y,2);%二阶导
y=dsolve(diff(y,3)-diff(y,2)==x,y(1)==9,dy(1)==7,d2y(2)==4);
y=simplify(y)%对结果进行化简

 想起了高数考试被手解常微分方程的恐惧😱。

更多内容参考 dsolve()文档。