1、.實驗01講評、參考答案講 評未交實驗報告的同學名單數學：01邊清水，09龔昱霏，14黃浦，34譚世韜信科：批改情況：不批改，同學們自己對照參考答案。附參考答案：數學建模實驗王平實驗01 建立數學模型（4學時）（第1章 建立數學模型）教材中給出原始數據，結合模型，得到結果。但如何求得結果這一過程沒有給出，實際上要用MATLAB軟件編寫程序來求得，這應該交給實驗課來完成。考慮到同學們剛學習MATLAB語言，編程能力不強，所以有關的程序給出來供同學們進行驗證。要求同學們要讀懂程序。1.（求解，編程）如何施救藥物中毒p1011人體胃腸道和血液系統中的藥量隨時間變化的規律（模型）：其中，x(t)為t時

2、刻胃腸道中的藥量，y(t)為t時刻血液系統中的藥量，t=0為服藥時刻。1.1（求解）模型求解p1011要求： 用MATLAB求解微分方程函數dsolve求解該微分方程（符號運算）。 用MATLAB的化簡函數simplify化簡所得結果。提示：dsolve和simplify的用法可用help查詢。建議在命令窗口中操作。 求解的語句及運行結果（比較11式(3)、(4)）：>> x,y= dsolve('Dx=-a*x','Dy=a*x-b*y','x(0)=1100','y(0)=0');>> disp(x,y

3、) 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)1.2（編程）結果分析p11已知=0.1386, =0.1155，將上題中得到x(t)和y(t)兩條曲線畫在同一個圖形窗口內（見11圖1）。提示：MATLAB命令：plot, fplot, hold on/off, grid on

4、/off, xlabel, ylabel, text。 編寫的程序和運行結果（比較11圖1）：程序1：用plotclc;clear;a=0.1386; b=0.1155;t=0:0.01:25;x=1100./exp(a*t);y=-(1100*a*(1./exp(a*t) -1./exp(b*t)/(a - b);plot(t,x,t,y);grid on;xlabel('itt /h'); ylabel('itx,ity /mg');text(2,1100/exp(a*2),' itx(itt)');text(3,-(1100*a*(1/ex

5、p(a*3) - 1/exp(b*3)/(a - b),' ity(itt)');程序2：用fplot和匿名函數clc;clear;a=0.1386; b=0.1155;fplot(t)1100/exp(a*t),-(1100*a*(1/exp(a*t) - 1/exp(b*t)/(a - b),0 25);grid on;xlabel('itt /h'); ylabel('itx,ity /mg');text(2,1100/exp(a*2),' itx(itt)');text(3,-(1100*a*(1/exp(a*3) - 1

6、/exp(b*3)/(a - b),' ity(itt)');2.（編程，驗證）商人們怎樣安全過河p89三名商人各帶一個隨從乘船渡河，一只小船只能容納二人，由他們自己劃行。隨從們密約，在河的任一岸，一旦隨從的人數比商人多，就殺人越貨。但是如何乘船的大權掌握在商人們手中。商人們怎樣才能安全渡河呢？模型構成決策： 每一步（此岸到彼岸或彼岸到此岸）船上的人員。要求：在安全的前提下（兩岸的隨從數不比商人多），經有限步使全體人員過河。xk第k次渡河前此岸的商人數yk第k次渡河前此岸的隨從數 xk , yk=0,1,2,3; k=1,2,¼過程的狀態sk=(xk , yk)允許狀

7、態集合S=(x, y)| x=0, y=0,1,2,3; x=3, y=0,1,2,3; x=y=1,2uk第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多步決策問題求dkÎD(k=1, 2, ¼, n), 使skÎS, 并按轉移律由 s1=(3,3) 到達sn+1=(0,0)。2.1（編程）求允許決策集合D和允許狀態集合SD是2行多列矩陣，每一列是一個決策。S是2行多列矩陣，每一列是一種狀態

8、。要求： 編寫一個命令文件的程序求D和S，并輸出。 S的第一列是3,3'，最后一列是0,0 '。 編寫的程序和運行結果：程序：clear; clc;%求允許決策集合D（2×n1，n1種決策）D=;for u=0:2 for v=0:2 if u+v=1|u+v=2 D=D,u;v; end endend%求允許狀態集合S（2×n2，n2種狀態）S=;for x=3:-1:0 for y=3:-1:0 if x=0|x=3|x=y S=S,x;y; end endend%首列狀態 ( 商人數, 仆從數)' = ( 3, 3 )'，末列為( 0,

9、 0 )'D, S運行結果：2.2（驗證）給出一個商人們安全過河的方案程序（輸入時，不必把注釋也輸入）：程序運行結果：要求： 讀懂以上程序。 將題2.1的程序放到上面的程序中（去掉最后多余的輸出語句），運行程序。對照答案，如不一致，請檢查程序。 給出運行的完整程序和運行結果（參考9圖1）：程序：function river()D,S1=DS(); S2=S1; SS=;%初始允許狀態集合S1,S2s1=S1(:,1); S1(:,1)=;%始發狀態s1=3;3，用過的狀態移出i=1; SS(:,i)=s1;%SS存放有效狀態轉移的路徑k=true;%k=true時，從此岸到彼岸；否則從

10、彼岸到此岸while isempty(S2)%目標狀態集合非空 s2,S2=fun(s1,S2,k,D); if k&&s2(1)=0&&s2(2)=0%有解 SS(:,i+1)=s2; disp(1:i+1;SS(:,1:i+1); return; end if s2(1)=-1%新狀態s2無效 if i=1%不能再回退 disp('無解！'); return; end i=i-1;%回退 s1=SS(:,i); else i=i+1; SS(:,i)=s2;%s1進入路徑 s1=s2; end S=S1;S1=S2;S2=S; k=k;end

11、disp('無解！');end function s2,S2=fun(s1,S2,k,D)%k=true時，從此岸到彼岸；否則從彼岸到此岸for u=D %決策遍歷 s2=s1+(-1)k*u; %狀態轉移s1->s2 for i=1:size(S2,2) if s2=S2(:,i)%是否新允許狀態 S2(:,i)=; return; end endends2=-1;-1;%從s1走不通end function D,S=DS()%求允許決策集合D（2×n1，n1種決策）D=;for u=0:2 for v=0:2 if u+v=1|u+v=2 D=D,u;v;

12、end endend%求允許狀態集合S（2×n2，n2種狀態）S=;for x=3:-1:0 for y=3:-1:0 if x=0|x=3|x=y S=S,x;y; end endend%首列狀態 ( 商人數, 仆從數)' = ( 3, 3 )'，末列為( 0, 0 )'end給出程序的運行結果：3.（求解）商人們怎樣安全過河（修改）p9對第2題的問題改動，用類似的方法求解。3.1 在第2題中修改商人數和隨從數有四名商人各帶一個隨從，其它同第2題。修改第2題中的程序求解。 修改的程序部分和完整程序的運行結果（安全過河的方案）：%求允許狀態集合S（2×

13、;n2，n2種狀態）S=;for x=4:-1:0 for y=4:-1:0 if x=0|x=4|x=y S=S,x;y; end endend%首列狀態 (商人數, 仆從數)' = ( 4, 4 )'，末列為( 0, 0 )'3.2 在3.1題中修改船容納的人數船能容納3人。修改3.1題中的程序求解。 修改的程序部分和完整程序的運行結果（安全過河的方案）：%求允許決策集合D（2×n1，n1種決策）D=;for u=0:3 for v=0:3 if u+v>=1 && u+v<=3 D=D,u;v; end endend4.（編程

14、）安全過河問題（人、貓、雞、米）p21習題5人帶著貓、雞、米過河，除需要人劃船之外，船至多能載貓、雞、米三者之一，而當人不在場時貓要吃雞、雞要吃米。試設計一個安全過河方案。模仿商人過河問題的程序編寫程序解決本問題。 程序：function rive()%求允許決策集合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;for a=23+22+21+1:-1:0 s=fix(a/23);mod(fix(a/22),2);mod(fix(a/2

15、),2);mod(a,2);%取二進制數字 if s(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)' S1=true(1,size(S,2);%S1為船在此岸時，當前未走過的狀態集合標記。 S2=S1;%S2（船在彼岸）類同。S1(1)=false;%狀態1做起始狀態SB=1;SE=;SS=;sn=-1;%狀態轉移SBSE，走過的“狀態轉移對”存入S

16、Swhile isempty(SB)%回溯法 if sn=-1, SM=S2; else SM=S1; end%確定目標狀態 for i=SB%狀態轉移SBSE for d=D sj=S(:,i)+sn*d; for j=1:size(S,2) if all(sj=S(:,j)&&SM(j)%有效的狀態并且沒走過 SS=SS,i;j;%將ij存入SS SE=SE,j;%到達的新狀態存入SE SM(j)=false;%標記已走過的狀態 if all(sj=S(:,end)%安全完成過河 show(SS,S); return;%輸出，程序運行結束 end break; end en

17、d end end if sn=-1, S2=SM; else S1=SM; end %更新 sn=-sn; SB=SE; SE=; %準備走下一步enddisp('無解！');end function show(SS,S)%輸出SD=SS(2,end); i=SS(1,end);for j=size(SS,2)-1:-1:1 if SS(2,j)=i SD=i,SD; i=SS(1,j); endendSD=1,SD;disp(1:length(SD);S(:,SD);end 程序的運行結果5.（選做，編程）商人們怎樣安全過河（自編）2.3中的程序有一定難度，所以給出了參考程

18、序，初學MATLAB的同學可能想不到用其中的一些語句。第2題給出的是求解商人們安全過河問題的一種算法程序，還有其它的算法程序，比如回溯法、遞歸算法等。要求： 按照自己的理解和想法，改寫2.3中的程序。 或用別的算法編寫解決商人們安全過河的程序。 程序和運行結果：解1（遞歸）程序function river()global mark D SHOW;%全局變量mark=2*ones(4);mark(1,:)=0; mark(4,:)=0;%mark(i,j)=0為允許狀態mark(2,2)=0; mark(3,3)=0;D=0 1; 0 2; 1 0; 1 1; 2 0;%允許決策集合mark(4

19、,4)=2; SHOW=;fun(3,3,-1);disp(3,3;SHOW);returnfunction f=fun(s,t)global mark D SHOW;%全局變量if s=0f=1; return;endfor k=1:5s1=s+t*D(k,:);i=s1(1)+1; j=s1(2)+1;if all(s1>=0) && all(s1<=3) && (mark(i,j)=0 | mark(i,j)=t)if mark(i,j)=0%標記狀態(i-1,j-1)的使用情況mark(i,j)=-t;elsemark(i,j)=2;end

20、if fun(s1,-t)SHOW=s1;SHOW; f=1; return; endendendif s=3 disp('無解');endf=0;運行結果3 33 13 23 0 %船停靠彼岸，先把所有隨從渡過彼岸3 11 12 20 2 %船停靠彼岸，再把所有商人渡過彼岸0 3 %船停靠此岸，此岸只有隨從，彼岸只有商人0 10 20 0 %船停靠彼岸，最后把所有隨從渡過彼岸三個階段：(3,3,船停此岸)(3,0,船停彼岸)(0,3,船停此岸)(0,0,船停彼岸)第一階段：所有隨從渡過彼岸（商人不動）；第二階段：使所有商人渡過彼岸，并將隨從渡回此岸；第三階段：所有隨從渡過彼岸

21、（商人不動）。解2（借用棧）思路定義狀態：此岸的商人數，隨從數，小船停靠在此岸或彼岸。從初始狀態開始，找到一個決策得到下一個狀態，再從新狀態開始，找到一個決策得到下一個狀態，.。在這個過程中，出現過的狀態將不可再用。若從某個狀態開始，所有決策都得不到新狀態，則返回上一狀態，并從下一個決策開始，試探是否能得到新的狀態。直到到達狀態(0,0)，或無解。算法1. 給出允許狀態集合；%用4階方陣mark標記，其中mark(i,j)=0表示(i-1,j-1)為允許狀態，mark(i,j)=2表示(i-1,j-1)為不可達狀態2. 給出允許決策集合；% D=0 1; 0 2; 1 0; 1 1; 2 03

22、. 初始狀態進棧；%包括：當前狀態，出發標記（-1/1為從此岸/彼岸出發），選擇第幾個決策4. 初始狀態標記；%mark(i,j)=-1/1表示從此岸/彼岸出發轉移到狀態(i-1,j-1)已走過，mark(i,j)=2表示上述兩種情況都出現或不可達狀態5. 循環直到棧空或到達狀態(0,0) 5.1 出棧，將出棧的狀態作為當前狀態，并選擇下一個決策；%原決策失敗，回溯 5.2 循環直到決策都用過或到達最終狀態 5.2.1 由當前狀態、出發標記和決策得到轉移到的新狀態； 5.2.2 若新狀態是允許狀態且沒走過，則執行5.2.2.1；否則選擇下一個決策 5.2.2.1 標記新狀態的訪問信息；%mar

23、k(i,j)=-1/1/2 5.2.2.2 當前狀態重新進棧（改變了決策）；%存儲決策過程 5.2.2.4 將新狀態作為當前狀態，置出發標記，選擇第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=0 1; 0 2; 1 0; 1 1; 2 0;%允許決策集合sp=1; st(sp,:)=3,3,-1,0; %進棧，(3,3)為當前狀態，-1表示小船停靠此岸，決策0（未做決策）mark(4,4)

24、=2;%標記為不可達狀態s0=3,3;while sp>0&& any(s0>0)%棧未空且s0中至少有一個元素>0s0=st(sp,1,2); ss=st(sp,3);%當前狀態，ss=-1/1小船停靠此岸/彼岸d=st(sp,4);%之前用過的決策，該決策失敗sp=sp-1;%出棧 d=d+1;%試探下一個決策 while d<6 && any(s0>0)s1=s0+ss*D(d,:);%下一個狀態，不包括小船停靠狀態i=s1(1)+1;j=s1(2)+1;%狀態轉換為下標if all(s1>=0) &&

25、all(s1<=3) && (mark(i,j)=0 | mark(i,j)=-ss)if mark(i,j)=0%標記狀態(i-1,j-1)的使用情況mark(i,j)=ss;elsemark(i,j)=2;endsp=sp+1;%進棧st(sp,:)=s0,ss,d;%當前狀態重新進棧，改變了決策s0=s1; ss=-ss; d=1;%下一個狀態作當前狀態，選擇第1個決策else%決策d失敗d=d+1;endendendif s0=0disp(1:sp+1',st(1:sp,1:2);s0);%有解elsedisp('無解'); end運行結果

26、 1 3 3 2 3 1 3 3 2 4 3 0 5 3 1 6 1 1 7 2 2 8 0 2 9 0 3 10 0 1 11 0 2 12 0 0解3clear; clc;%求允許狀態集合SS=;for x=0:3 for y=0:3 if (x=0&&(y>=0&&y<=3). |(x=3&&(y>=0&&y<=3). |(x=1&&y=1|x=2&&y=2) S=S,x;y; end endendSA=S;%從此岸(A)到彼岸(B)渡河前此岸的允許狀態集合SB=S;%

27、從彼岸(B)到此岸(A)渡河前此岸的允許狀態集合%求允許決策集合DD=;for u=0:2 for v=0:2 if u+v=1|u+v=2 D=D,u;v; end endend k=1;%第1次渡河s=3;3;%第1次渡河前的狀態Snum(1)=1;%第1次渡河前的狀態總數SA(:,end)=; SB=SA;%去掉狀態(3,3)while Snum(k)>0 Snum(k+1)=0;%第k+1次渡河前的狀態總數初始化 k1=size(s,2)-Snum(k)+1;k2=size(s,2); for i=k1:k2 xk=s(1,i); yk=s(2,i); for d=D xk1=x

28、k+(-1)k*d(1); yk1=yk+(-1)k*d(2); if (-1)k=1 for j=1:size(SB,2) if xk1=SB(1,j)&&yk1=SB(2,j) s=s,SB(:,j); SB(:,j)=; Snum(k+1)=Snum(k+1)+1; break; end end else for j=1:size(SA,2) if xk1=SA(1,j)&&yk1=SA(2,j) s=s,SA(:,j); SA(:,j)=; Snum(k+1)=Snum(k+1)+1; break; end end end end end k=k+1;e

29、ndSnums(:,1:Snum(1)m=Snum(1);for k=1:11 disp(s(:,m+1:m+Snum(k+1); m=m+Snum(k+1);end解4i=1;S=;D=;S(:,i)=3;3;%開始S(:,1)=3,3'while S(1,i)=0|S(2,i)=0%當s(:,i)=0,0'完成 D(:,i)=round(2*rand(2,1);%隨機取0,1,2給決策D(:,i) d=D(:,i); if d(1)+d(2)=1&&d(1)+d(2)=2%成立則不是允許決策 continue;%重新選取決策 end S(:,i+1)=S(:

30、,i)+(-1)i*D(:,i);%得到下一個狀態S(:,i+1) s=S(:,i+1); if (s(1)=1&&s(2)=1)|(s(1)=2&&s(2)=2)|(s(1)>3|s(1)<0)|(s(2)>3|s(2)<0)%成立則不是允許狀態 continue;%重新選取決策 end i=i+1; if i=105%等于此循環次數判無解 disp('無解');return; endendfor i=1:size(S,2)%添加狀態的第3行元素，-1為船在此岸，1時在彼岸 SS(:,i)=S(:,i);(-1)i;en

31、d%區分過程中的不同狀態for i=1:size(SS,2)-1%去掉狀態轉移過程中多余的狀態 a=0; for j=i+1:size(SS,2)%找出與第i狀態相同的最后一個狀態 if any(SS(:,i)-SS(:,j)=0 a=j; end end if i<a%刪除第i到a-1列，多余的狀態轉移 SS(:,i:a-1)=; endenddisp(1:size(SS,2);SS)%該程序運行效率含隨機性，可能誤判為無解！解5syms k;%定義一個關于k的，儲存各個狀態之間轉化的符號矩陣as=3,2,1,0;ss=;%能夠存在的狀態集合for i=1:length(s)%計算ss

32、 for j=1:length(s) if i=length(s)|i=1|s(i)=s(j) sss=s(i),s(j); ss=ss;sss; end endendignore,p=sort(sum(ss,2),'descend');%按總數排序,便于計算狀態矩陣ss=ss(p,:);a=sym(zeros(length(ss),length(ss);for i=1:length(ss) for j=i+1:length(ss) m=abs(sum(ss(i,:)-sum(ss(j,:); if m=1|m=2%兩個狀態的總人數之差只能是1或者2 if ss(i,:)>

33、;=ss(j,:)%保證總人數多的狀態每個分量都要大 a(i,j)=2+(-1)k; a(j,i)=(-1)k; end end endenda(1,2)=0;%當第一次由(3,3)變到(3,2)時，下一次必為(3,2) 變到(3,3)，避免無用功k=1;dd=1;s0=ss(1,:);%初始狀態S=ss(1,:);%儲存總的狀態轉移方案D=;%每一步的決策集合while any(s0)&k<=100 b=eval(a);%確定此時的a矩陣 for i=1:length(ss)%判斷s0處于s中的第幾行 if s0(1,:)=ss(i,:) break end end m=fin

34、d(b(i,:)=1); %確定此時的方案集合 if k>1&length(m)>=2 m(find(m=dd(k-1)=;%避免走回頭路 end ignore,p=sort(rand(1,length(m);%在剩下的狀態中隨機找一個可行狀態 s0=ss(m(p(1),:); k=k+1; dd=dd;m(p(1);S=S;s0;D=D;abs(S(k,:)-S(k-1,:);endSD解6clear; clc;%求允許決策集合D（2×n1，n1種決策）D=;for u=0:2 for v=0:2 if u+v=1|u+v=2 D=D,u;v; end ende

35、nd%求允許狀態集合S（2×n2，n2種狀態）S=;for x=3:-1:0 for y=3:-1:0 if x=0|x=3|x=y S=S,x;y; end endend%首列狀態 ( 商人數, 仆從數)' = ( 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,-

36、1'（起點），SS(:,end)=0,0,1'（終點）%狀態轉移矩陣A，A(i,j)=1表示存在決策使狀態i轉到j，其它為0A=zeros(SSnum);for i=1:SSnum for j=1:SSnum for d=D%順序取D的每一列給d s=SS(1:2,i)+SS(3,i)*d;-SS(3,i); if all(s=SS(:,j)%所有元素不為0時為真 A(i,j)=1; break; end end endend%SK為多行SSnum列的矩陣，第i列對應SS的第i列的狀態。%第i行表示第i次渡河前有那些狀態，1表示有對應下標列號的狀態，0則無。%第i+1行是第i行

37、狀態轉移過來的狀態%當第SSnum狀態出現時停止，表示有解。%或者新行的元素全0時停止，表示無解。%新行不允許出現之前出現過的狀態。s=1,zeros(1,SSnum-1);%從狀態3,3,-1'開始SK=s;while any(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%多步決策，直到找到目標。有難度！if any(s)=0 %新行全0disp('無解！'); return;end%最后一行的第SSnum狀態(j)開始往回找轉換到它的狀態(i)，直到狀態1。%得到一個狀態號的順序sk（行向量），為一個解。sk=zeros(1,size(SK,1);sk(1)=1; sk(1,end)=SSnum; j=SSnum;for k=size(SK,1)-1:-1:2 for i=find(SK(k,:) if A(i,j)=1 sk(k)=i; j=i; break; end endenddisp(1:length(sk);S