数学模型程序代码Matlab姜启源第一章建立数学模型.docx
《数学模型程序代码Matlab姜启源第一章建立数学模型.docx》由会员分享，可在线阅读，更多相关《数学模型程序代码Matlab姜启源第一章建立数学模型.docx（57页珍藏版）》请在冰豆网上搜索。
数学模型程序代码Matlab姜启源第一章建立数学模型
第1章建立数学模型
教材中给出原始数据，结合模型，得到结果。
但如何求得结果这一过程没有给出，实际上要用MATLAB软件编写程序来求得，这应该交给实验课来完成。
考虑到上学期同学们刚学习MATLAB语言，编程能力不强，所以有关的程序给出来供同学们进行验证。
要求同学们要读懂程序。
1.（求解，编程）如何施救药物中毒p10~11
人体胃肠道和血液系统中的药量随时间变化的规律（模型）：
其中，x（t）为t时刻胃肠道中的药量，y（t）为t时刻血液系统中的药量，t=0为服药时刻。
（求解）模型求解p10~11
要求：
①用MATLAB求解微分方程函数dsolve求解该微分方程（符号运算）。
②用MATLAB的化简函数simplify化简所得结果。
③结果与教材P11上的内容比较。
提示：
dsolve和simplify的用法可用help查询。
建议在命令窗口中操作。
★求解的语句及运行结果：
>>[x,y]=dsolve（'Dx=-a*x','Dy=a*x-b*y','x（0）=1100','y（0）=0'）;
>>disp（[x,y]）
[1100*exp（-a*t）,exp（-a*t）*exp（-b*t）*（（1100*a*exp（a*t））/（a-b）-（1100*a*exp（b*t））/（a-b））]
>>disp（simplify（[x,y]））;
[1100*exp（-a*t）,（1100*a*exp（-t*（a+b））*（exp（a*t）-exp（b*t）））/（a-b）]
（编程）结果分析p11
已知λ=,μ=，将上题中得到x（t）和y（t）两条曲线画在同一个图形窗口内。
参考图形如下。
提示：
MATLAB命令plot,fplot,holdon/off,gridon/off,xlabel,ylabel,text。
★编写的程序和运行结果：
程序1：
用plot
clc;clear;
a=;b=;
t=0:
:
25;
x=1100./exp（a*t）;
y=-（1100*a*（1./exp（a*t）-1./exp（b*t）））/（a-b）;
plot（t,x,t,y）;
gridon;
xlabel（'{\itt}/h'）;ylabel（'{\itx},{\ity}/mg'）;
text（2,1100/exp（a*2）,'{\itx}（{\itt}）'）;
text（3,-（1100*a*（1/exp（a*3）-1/exp（b*3）））/（a-b）,'{\ity}（{\itt}）'）;
程序2：
用fplot和匿名函数
clc;clear;
a=;b=;
fplot（@（t）[1100/exp（a*t）,-（1100*a*（1/exp（a*t）-1/exp（b*t）））/（a-b）],[025]）;
gridon;
xlabel（'{\itt}/h'）;ylabel（'{\itx},{\ity}/mg'）;
text（2,1100/exp（a*2）,'{\itx}（{\itt}）'）;
text（3,-（1100*a*（1/exp（a*3）-1/exp（b*3）））/（a-b）,'{\ity}（{\itt}）'）;
2.（编程，验证）商人们怎样安全过河p8~9
三名商人各带一个随从乘船渡河，一只小船只能容纳二人，由他们自己划行。
随从们密约，在河的任一岸，一旦随从的人数比商人多，就杀人越货。
但是如何乘船的大权掌握在商人们手中。
商人们怎样才能安全渡河呢
[模型构成]
决策：
每一步（此岸到彼岸或彼岸到此岸）船上的人员。
要求：
在安全的前提下（两岸的随从数不比商人多），经有限步使全体人员过河。
xk第k次渡河前此岸的商人数
yk第k次渡河前此岸的随从数xk,yk=0,1,2,3;k=1,2,
过程的状态sk=（xk,yk）
允许状态集合S={（x,y）x=0,y=0,1,2,3;x=3,y=0,1,2,3;x=y=1,2}
uk第k次渡船上的商人数
vk第k次渡船上的随从数uk,vk=0,1,2;k=1,2,
决策dk=（uk,vk）
允许决策集合D={（u,v）u+v=1,2}
状态转移律sk+1=sk+（-1）kdk
[多步决策问题]
求dkD（k=1,2,,n）,使skS,并按转移律由s1=（3,3）到达sn+1=（0,0）。
（编程）求允许决策集合D和允许状态集合S
D是2行多列矩阵，每一列是一个决策。
S是2行多列矩阵，每一列是一种状态。
要求：
①编写程序求D和S，并输出。
②S的第一列是[3,3]'，最后一列是[0,0]'。
★编写的程序和运行结果：
程序：
clear;clc;
%求允许决策集合D（2×n1，n1种决策）
D=[];
foru=0:
2
forv=0:
2
ifu+v==1||u+v==2
D=[D,[u;v]];
end
end
end
%求允许状态集合S（2×n2，n2种状态）
S=[];
forx=3:
-1:
0
fory=3:
-1:
0
ifx==0||x==3||x==y
S=[S,[x;y]];
end
end
end%首列状态（商人数,仆从数）'=（3,3）'，末列为（0,0）'
D,S
运行结果：
（验证）求动态允许状态集合SS和状态转移矩阵A
上面允许状态集合S没有指明当时船是在此岸还是在彼岸，应该将S中的每一种状态再分为两种状态，需增加一个元素（值为-1或1）放在第三行。
定义动态允许状态集合
SS={（x,y,z）'x=0,y=0,1,2,3;x=3,y=0,1,2,3;x=y=1,2;z=-1,1}
（x,y,-1），表示从此岸渡河前此岸的允许状态（x,y）。
（x,y,1），表示从彼岸渡河前此岸的允许状态（x,y）。
SS是三行多列矩阵，每一列表示一种状态，列下标为其编号。
定义状态转移矩阵A，其中，A（i,j）=1表示D中存在决策使状态i转到j，否则A（i,j）=0。
程序如下（输入时，不必把注释也输入）：
程序运行结果（参考）：
要求：
①将程序接在上题的程序之后（去掉最后多余的输出语句），程序最后给出显示SS和A的语句。
②运行程序，输出SS和A。
对照参考答案，如数值不一致，请检查程序。
MATLAB函数all的用法见提示。
★运行的完整程序和运行结果：
程序：
clear;clc;
%求允许决策集合D（2×n1，n1种决策）
D=[];
foru=0:
2
forv=0:
2
ifu+v==1||u+v==2
D=[D,[u;v]];
end
end
end
%求允许状态集合S（2×n2，n2种状态）
S=[];
forx=3:
-1:
0
fory=3:
-1:
0
ifx==0||x==3||x==y
S=[S,[x;y]];
end
end
end%首列状态（商人数,仆从数）'=（3,3）'，末列为（0,0）'
%动态允许状态集合SS（3×n3，n3种状态）
%-1，从此岸渡河前此岸的允许状态
%1，从彼岸渡河前此岸的允许状态
SS=[[S;-ones（1,size（S,2））],[S;ones（1,size（S,2））]];
SSnum=size（SS,2）;%状态总数，SS中的列下标对应状态的编号
%SS（:
1）=[3,3,-1]'（起点），SS（:
end）=[0,0,1]'（终点）
%状态转移矩阵A，A（i,j）=1表示存在决策使状态i转到j，其它为0
A=zeros（SSnum）;
fori=1:
SSnum
forj=1:
SSnum
ford=D%顺序取D的每一列给d
s=[SS（1:
2,i）+SS（3,i）*d;-SS（3,i）];
ifall（s==SS（:
j））%所有元素不为0时为真
A（i,j）=1;break;
end
end
end
end
SS,A
运行结果：
（验证）给出一个商人们安全过河的方案
程序如下（输入时，不必把注释也输入）：
程序运行结果：
要求：
①读懂以上程序，掌握all和any函数的应用。
②将程序接在之前的程序后（去掉最后多余的输出语句），运行程序。
对照答案，如不一致，请检查程序。
MATLAB函数any的用法见提示。
★给出运行的完整程序和运行结果：
程序：
clear;clc;
%求允许决策集合D（2×n1，n1种决策）
D=[];
foru=0:
2
forv=0:
2
ifu+v==1||u+v==2
D=[D,[u;v]];
end
end
end
%求允许状态集合S（2×n2，n2种状态）
S=[];
forx=3:
-1:
0
fory=3:
-1:
0
ifx==0||x==3||x==y
S=[S,[x;y]];
end
end
end%首列状态（商人数,仆从数）'=（3,3）'，末列为（0,0）'
%动态允许状态集合SS（3×n3，n3种状态）
%-1，从此岸渡河前此岸的允许状态
%1，从彼岸渡河前此岸的允许状态
SS=[[S;-ones（1,size（S,2））],[S;ones（1,size（S,2））]];
SSnum=size（SS,2）;%状态总数，SS中的列下标对应状态的编号
%SS（:
1）=[3,3,-1]'（起点），SS（:
end）=[0,0,1]'（终点）
%状态转移矩阵A，A（i,j）=1表示存在决策使状态i转到j，其它为0
A=zeros（SSnum）;
fori=1:
SSnum
forj=1:
SSnum
ford=D%顺序取D的每一列给d
s=[SS（1:
2,i）+SS（3,i）*d;-SS（3,i）];
ifall（s==SS（:
j））%所有元素不为0时为真
A（i,j）=1;break;
end
end
end
end
%SK为多行SSnum列的矩阵，第i列对应SS的第i列的状态。
%第i行表示第i次渡河前有那些状态，1表示有对应下标列号的状态，0则无。
%第i+1行是第i行状态转移过来的状态
%当第SSnum状态出现时停止，表示有解。
%或者新行的元素全0时停止，表示无解。
%新行不允许出现之前出现过的状态。
s=[1,zeros（1,SSnum-1）];%从状态[3,3,-1]'开始
SK=s;
whileany（s）&&s（SSnum）~=1
k=find（s）;%求最后一行为1的列下标
s=any（A（k,:
）,1）;%得到新行，元素值为1或0
s=s&~any（SK,1）;%去掉新行中之前出现的状态
SK=[SK;s];%新行加入，作为SK的末行
end%多步决策，直到找到目标。
有难度！
ifany（s）==0%新行全0
disp（'无解！
'）;return;
end
%最后一行的第SSnum状态（j）开始往回找转换到它的状态（i），直到状态1。
%得到一个状态号的顺序sk（行向量），为一个解。
sk=zeros（1,size（SK,1））;
sk
（1）=1;sk（1,end）=SSnum;j=SSnum;
fork=size（SK,1）-1:
-1:
2
fori=find（SK（k,:
））
ifA（i,j）==1
sk（k）=i;j=i;break;
end
end
end
disp（[1:
length（sk）;SS（1:
2,sk）]）;%输出结果
给出程序的运行结果：
3.（求解）商人们怎样安全过河（修改）p9
对第2题的问题改动，用类似的方法求解。
在第2题中修改商人数和随从数
有四名商人各带一个随从，其它同第2题。
修改第2题中的程序求解。
★修改的程序部分和完整程序的运行结果（安全过河的方案）：
%求允许状态集合S（2×n2，n2种状态）
S=[];
forx=4:
-1:
0
fory=4:
-1:
0
ifx==0||x==4||x==y
S=[S,[x;y]];
end
end
end%首列状态（商人数,仆从数）'=（4,4）'，末列为（0,0）'
在题中修改船容纳的人数
船能容纳3人。
修改题中的程序求解。
★修改的程序部分和完整程序的运行结果（安全过河的方案）：
%求允许决策集合D（2×n1，n1种决策）
D=[];
foru=0:
3
forv=0:
3
ifu+v>=1&&u+v<=3
D=[D,[u;v]];
end
end
end
4.（编程）安全过河问题（人、猫、鸡、米）p21习题5
人带着猫、鸡、米过河，除需要人划船之外，船至多能载猫、鸡、米三者之一，而当人不在场时猫要吃鸡、鸡要吃米。
试设计一个安全过河方案。
模仿商人过河问题的程序编写程序解决本问题。
★程序：
clear;clc;
%求允许决策集合D（2×n1，n1种决策）
D=[[1,0,0,0]',[1,1,0,0]',[1,0,1,0]',[1,0,0,1]'];
%求允许状态集合S（2×n2，n2种状态）
S=[];j=1;
fora=2^3+2^2+2^1+1:
-1:
0
s=[fix（a/2^3）;mod（fix（a/2^2）,2）;mod（fix（a/2）,2）;mod（a,2）];%取二进制数字
ifs
（1）&&（s
（2）&&s（4）||s（3））||~s
（1）&&（~s
（2）&&~s（4）||~s（3））
S（:
j）=s;j=j+1;
end%s=[人，猫，鸡，米]
end%首列状态=（1,1,1,1）'，末列为（0,0,0,0）'
%动态允许状态集合SS（3×n3，n3种状态）
%-1，从此岸渡河前此岸的允许状态
%1，从彼岸渡河前此岸的允许状态
SS=[[S;-ones（1,size（S,2））],[S;ones（1,size（S,2））]];
SSnum=size（SS,2）;%状态总数，SS中的列下标对应状态的编号
%SS（:
1）=[1,1,1,1,-1]'（起点）,SS（:
end）=[0,0,0,0,1]'（终点）
%状态转移矩阵A，A（i,j）=1表示存在决策使状态i转到j，其它为0
A=zeros（SSnum）;
fori=1:
SSnum
forj=1:
SSnum
ford=D%顺序取D的每一列给d
s=[SS（1:
4,i）+SS（5,i）*d;-SS（5,i）];
ifall（s==SS（:
j））%所有元素不为0时为真
A（i,j）=1;break;
end
end
end
end
%SK为多行SSnum列的矩阵，第i列对应SS的第i列的状态。
%第i行表示第i次渡河前有那些状态，1表示有对应下标列号的状态，0则无。
%第i+1行是第i行状态转移过来的状态
%当第SSnum状态出现时停止，表示有解。
%或者新行的元素全0时停止，表示无解。
%新行不允许出现之前出现过的状态。
s=[1,zeros（1,SSnum-1）];%从状态[1,1,1,1,-1]'开始
SK=s;
whileany（s）&&all（SSnum-find（s））
k=find（s）;%求最后一行为1的列下标
s=any（A（k,:
）,1）;%得到新行，元素值为1或0
s=s&~any（SK,1）;%去掉新行中之前出现的状态
SK=[SK;s];%新行加入，作为SK的末行
end
if~any（s）%新行全0
disp（'无解！
'）;return;
end
%最后一行的第SSnum状态（j）开始往回找转换到它的状态（i），直到状态1。
%得到一个状态号的顺序sk（行向量），为一个解。
sk=zeros（1,size（SK,1））;
sk
（1）=1;sk（1,end）=SSnum;j=SSnum;
fork=size（SK,1）-1:
-1:
2
fori=find（SK（k,:
））
ifA（i,j）==1
sk（k）=i;j=i;break;
end
end
end
disp（[1:
length（sk）;SS（1:
4,sk）]）;%输出结果
★程序的运行结果
5.（选做，编程）商人们怎样安全过河（自编）
中的程序有一定难度，所以给出了参考程序，初学MATLAB的同学可能想不到用其中的一些语句。
第2题给出的是求解商人们安全过河问题的一种算法程序，还有其它的算法程序，比如回溯法、递归算法等。
要求：
①按照自己的理解和想法，改写中的程序。
②或用别的算法编写解决商人们安全过河的程序。
★程序和运行结果：
解1（递归）
程序
functionriver（）
globalmarkDSHOW;%全局变量
mark=2*ones（4）;
mark（1,:
）=0;mark（4,:
）=0;%mark（i,j）=0为允许状态
mark（2,2）=0;mark（3,3）=0;
D=[01;02;10;11;20];%允许决策集合
mark（4,4）=2;SHOW=[];
fun（[3,3],-1）;
disp（[[3,3];SHOW]）;
return
functionf=fun（s,t）
globalmarkDSHOW;%全局变量
ifs==0
f=1;return;
end
fork=1:
5
s1=s+t*D（k,:
）;
i=s1
（1）+1;j=s1
（2）+1;
ifall（s1>=0）&&all（s1<=3）&&（mark（i,j）==0||mark（i,j）==t）
ifmark（i,j）==0%标记状态（i-1,j-1）的使用情况
mark（i,j）=-t;
else
mark（i,j）=2;
end
iffun（s1,-t）
SHOW=[s1;SHOW];f=1;
return;
end
end
end
ifs==3
disp（'无解'）;
end
f=0;
运行结果
33
31
32
30%船停靠彼岸，先把所有随从渡过彼岸
31
11
22
02%船停靠彼岸，再把所有商人渡过彼岸
03%船停靠此岸，此岸只有随从，彼岸只有商人
01
02
00%船停靠彼岸，最后把所有随从渡过彼岸
三个阶段：
（3,3,船停此岸）→（3,0,船停彼岸）→（0,3,船停此岸）→（0,0,船停彼岸）
第一阶段：
所有随从渡过彼岸（商人不动）；
第二阶段：
使所有商人渡过彼岸，并将随从渡回此岸；
第三阶段：
所有随从渡过彼岸（商人不动）。
解2（借用栈）
思路
定义状态：
此岸的商人数，随从数，小船停靠在此岸或彼岸。
从初始状态开始，找到一个决策得到下一个状态，再从新状态开始，找到一个决策得到下一个状态，...。
在这个过程中，出现过的状态将不可再用。
若从某个状态开始，所有决策都得不到新状态，则返回上一状态，并从下一个决策开始，试探是否能得到新的状态。
直到到达状态（0,0），或无解。
算法
1.给出允许状态集合；%用4阶方阵mark标记，其中mark（i,j）=0表示（i-1,j-1）为允许状态，mark（i,j）=2表示（i-1,j-1）为不可达状态
2.给出允许决策集合；%D=[01;02;10;11;20]
3.初始状态进栈；%包括：
当前状态，出发标记（-1/1为从此岸/彼岸出发），选择第几个决策
4.初始状态标记；%mark（i,j）=-1/1表示从此岸/彼岸出发转移到状态（i-1,j-1）已走过，mark（i,j）=2表示上述两种情况都出现或不可达状态
5.循环直到栈空或到达状态（0,0）
出栈，将出栈的状态作为当前状态，并选择下一个决策；%原决策失败，回溯
循环直到决策都用过或到达最终状态
由当前状态、出发标记和决策得到转移到的新状态；
若新状态是允许状态且没走过，则执行；否则选择下一个决策
标记新状态的访问信息；%mark（i,j）=-1/1/2
当前状态重新进栈（改变了决策）；%存储决策过程
将新状态作为当前状态，置出发标记，选择第1个决策；
6.若有解，则输出决策过程；否则输出“无解”；
程序
clear;clc;
mark=2*ones（4）;
mark（1,:
）=0;mark（4,:
）=0;%mark（i,j）=0为允许状态
mark（2,2）=0;mark（3,3）=0;
D=[01;02;10;11;20];%允许决策集合
sp=1;st（sp,:
）=[3,3,-1,0];%进栈，（3,3）为当前状态，-1表示小船停靠此岸，决策0（未做决策）
mark（4,4）=2;%标记为不可达状态
s0=[3,3];
whilesp>0&&any（s0>0）%栈未空且s0中至少有一个元素>0
s0=st（sp,[1,2]）;ss=st（sp,3）;%当前状态，ss=-1/1小船停靠此岸/彼岸
d=st（sp,4）;%之前用过的决策，该决策失败
sp=sp-1;%出栈
d=d+1;%试探下一个决策
whiled<6&&any（s0>0）
s1=s0+ss*D（d,:
）;%下一个状态，不包括小船停靠状态
i=s1
（1）+1;j=s1
（2）+1;%状态转换为下标
ifall（s1>=0）&&all（s1<=3）&&（mark（i,j）==0||mark（i,j）==-ss）
ifmark（i,j）==0%标记状态（i-1,j-1）的使用情况
mark（i,j）=ss;
else
mark（i,j）=2;
end
sp=sp+1;%进栈
st（sp,:
）=[s0,ss,d];%当前状态重新进栈，改变了决策
s0=s1;ss=-ss;d=1;%下一个状态作当前状态，选择第1个决策
else%决策d失败
d=d+1;
end
end
end
ifs0==0
disp（[[1:
sp+1]',[st（1:
sp,1:
2）;s0]]）;%有解
else
disp（'无解'）;
end
运行结果
133
231
332
430
531
611
722
802
903
1001
1102
升级会员 