前言
主要是通过matlab求解微分方程的实例,先来看一个小小问题:
\[x^2+y+(x-2y)y^{'}=0\]
可以通过matlab
的dsolve
函数求得通解,两行代码就能够解决问题,关于dsolve
的具体用法,参考官方文档的介绍
clc,clear
syms x %定义符号变量
y=dsolve('x^2+y+(x-2*y)*Dy==0','x')
%运行结果
% y =
% x/2 + ((4*x^3)/3 + x^2 + C1)^(1/2)/2
% x/2 - ((4*x^3)/3 + x^2 + C1)^(1/2)/2
% C1是常数
问题一
求解常微分方程组的通解
代码:
clc,clear
syms f(x) g(x)
[f1,g1]=dsolve('D2f+3*g==sin(x)','Dg+Df==cos(x)');
f1=simplify(f1) %化简
g1=simplify(g1)
% 运行结果
% f1 =
% C2 - sin(x)/3 + t*cos(x) + (3^(1/2)*C3*exp(3^(1/2)*t))/3 - (3^(1/2)*C4*exp(-3^(1/2)*t))/3
% g1 =
%
% sin(x)/3 - (3^(1/2)*C3*exp(3^(1/2)*t))/3 + (3^(1/2)*C4*exp(-3^(1/2)*t))/3
初值条件:
\[f^{'}(2)=0,f(3)=3,g(5)=1\]
求特解:
clc,clear
syms f(x) g(x)
[f2,g2]=dsolve('D2f+3*g==sin(x)','Dg+Df==cos(x)','Df(2)==0','f(3)==3','g(5)==1');
f2=simplify(f2)
g2=simplify(g2)
问题二
求解非齐次线性微分方程组:
\[X'=AX+f(t),X(t_0)=X_0\]
的解为:
\[X(t)=e^{A(t-t_0)}X_0+\int_{t_0}^{t} e^{A(t-s)}f(s)ds\]
clc,clear
syms x(t) y(t) z(t)
X=[x;y;z];
A=[1,0,0;2,1,-2;3,2,1];
B=[0;0;exp(t)*cos(2*t)];
X0=[0;1;1]; %初值条件
X=dsolve(diff(X)==A*X+B,X(0)==X0);
X=simplify([X.x;X.y;X.z])
% 运行结果
% X =
% 0
% -(exp(t)*(2*sin(2*t) - 2*cos(2*t) + t*sin(2*t)))/2
% (exp(t)*(4*cos(2*t) + 5*sin(2*t) + 2*t*cos(2*t)))/4
问题四:层次分析法(AHP)
主要用于解决评价类的问题,比如今天吃什么,明天去哪里旅游这类的。具体的开始通过案例来学习了解,更专业的介绍很多数模的书上都有,就不再此过多说明。
案例:Alice思考今天中午去哪里吃饭,有四家店A(饺子馆)、B(川菜馆)、C(KFC)、D(日料)可供选择。Alice考虑的重点因素有味道、环境、价格、服务。那么应该选择那家店享用今日的午餐呢?
1、建立层次结构模型
2、构造判断矩阵
3、层次但排序及其一致性检验
4、层次总排序及其一致性检验
建立层次结构模型
构造判断矩阵
准则层
B1 味道 | 1 | 3 | 2 | 1/2 |
B2 价格 | 1/3 | 1 | 1/2 | 1/4 |
B3 环境 | 1/2 | 2 | 1 | 1/3 |
B4 服务 | 2 | 4 | 3 | 1 |
方案层
A | 1 | 1/2 | 2 | 1/3 | A | 1 | 5 | 3 | 6 |
B | 2 | 1 | 5 | 3 | B | 1/5 | 1 | 1/2 | 4 |
C | 1/2 | 1/5 | 1 | 1/4 | C | 1/3 | 2 | 1 | 5 |
D | 3 | 1/3 | 4 | 1 | D | 1/6 | 1/4 | 1/5 | 1 |
B3 环境 | A | B | C | D | B4 服务 | A | B | C | D |
A | 1 | 1/2 | 1/2 | 1/4 | A | 1 | 1/4 | 1/3 | 1/6 |
B | 2 | 1 | 1/3 | 1/3 | B | 4 | 1 | 3 | 1/3 |
C | 2 | 3 | 1 | 1/2 | C | 3 | 1/3 | 1 | 1/5 |
D | 4 | 3 | 2 | 1 | D | 6 | 3 | 5 | 1 |
Matlab求解
clc,clear
% 读取txt中原始数据
data = fopen('1.txt','r');
% 准则数量
n1=4;
%方案数量
n2=4;
a=[];
for i=1:n1
%str2num:把字符串转换成数
%fgetl:用于从文件中读取一行数据,并丢弃行末的换行符
temp=str2num(fgetl(data));
% 准则层判断矩阵
a=[a;temp];
end
for i=1:n1
%ins2str:整形转换成字符串
% eval:将括号内的字符串视为语句并运行
str1=char(['b',int2str(i),'=[];']);
str2=char(['b',int2str(i),'=[b',int2str(i),';temp];']);
eval(str1);
for j =1:n2
temp=str2num(fgetl(data));
%读方案层判断矩阵
eval(str2);
end
end
ri=[0,0,0.58,0.90,1.12,1.24,1.32,1.41,1.45];%一致性指标
% eig:计算矩阵的特征值和特征向量
% 返回矩阵A特征值的对角矩阵 y 和矩阵 x
[x,y]=eig(a);
% 返回 y 的主对角线元素的列向量的最大值
lamda=max(diag(y));
% find:查找非零元素的索引和值
num=find(diag(y)==lamda);
w0=x(:,num)/sum(x(:,num));
% 准则层一致性检验
cr0=(lamda-n1)/(n1-1)/ri(n1)
if cr0<0.10
disp('通过一致性检验!');
end
for i=1:n1
[x,y]=eig(eval(char(['b',int2str(i)])));
lamda=max(diag(y));
num=find(diag(y)==lamda);
w1(:,i)=x(:,num)/sum(x(:,num));
% 方案层一致性检验
cr1=(lamda-n2)/(n2-1)/ri(n2);
end
cr1
if cr1<0.10
disp('通过一致性检验!');
end
% ts:准则层权重 cr:方案总权重
ts=w1*w0,cr=cr1*w0
运行结果:
cr0 = 0.0115 < 0.1
cr1 = 0.0545 < 0.1
ts =
0.1423
0.2909
0.1468
0.4200
cr =
0.0151
0.0052
0.0087
0.0255
所以Alice应该选择去吃日料