宏觀課件matlab課堂程序例題

1、本節(jié)課的主要目的在于讓了解一些的基本知識。因為大分部同學(xué)都沒有怎么接觸過這個，所以前半節(jié)我會講一些簡單的操作，大家跟著做下，基本上就能明白。而后半節(jié)課內(nèi)容會難一點，但我也會盡力讓大家明白的有趣之處。方面：的安裝與運行。沒有的同學(xué)可以下課來找我，然后我默認大家都會安裝這個，并且能順利打開它。打開之后，首先看見的界面叫做“命令窗口”，光標閃的地方就是你編輯命令的地方，在這個窗口每一句輸入的指令都會立即執(zhí)行。（嘗試下 clc）而點擊界面左上角第一個圖標可以新建一個程序，界面像一張空白頁，如下圖所示，編入這里令不會立即執(zhí)行，并且可以讓你修改，所以這里才是程序編輯的地方。編程完成之后，你需要保存該文件，

2、操作是 FileSave as選擇保存路徑取個好名字哦保存。這樣一來你就生成了一個后綴名為.m 的文件，這類文件就被叫做M 文件。之后，如何重新打開它們，例如打開 XX.m。操作是 FileOpen選擇路徑XX.m(XX:剛才你輸入的那個文件名)。最后，在子程序的運行操作中點擊圖標或者快捷鍵”F5”便可運行程序，結(jié)果將會出現(xiàn)在之前令窗口中。的函數(shù)和命令有很多，大家想深入了解，可以有以下方法：如果你知道指令，但不清楚用法，可以使用自帶的 help 命令 如：help global如果你不知道指令，那可以進行模糊搜索，在命令窗口輸入相關(guān)詞按 Tab 鍵就能出現(xiàn)所有相關(guān)令；3.當(dāng)然你也可以“加相關(guān)指

3、令” 如：“clf”。學(xué)= =# 如果你連指令是什么都不是很清楚的話，就先把這下吧。（講義中，表示在命令窗口輸入令，而%表示的是注釋，*后面的表示有難度可以跳過）的特點高效方便的矩陣數(shù)組運算靈活簡單的繪圖功能庫函數(shù)豐富，避免了繁雜的子程序編程任務(wù)*眾多面向領(lǐng)域應(yīng)用的工具箱主要介紹的是前面的 3 個特點，而最后的工具箱主要面對較高水平的應(yīng)用，所以本課是不會涉及到的。編輯程序：其實編程就是人機交流，計算機有自己能識別的語言，所以你得用這些計算機語言來告訴它是什么，怎么做。 這里，“是什么”即賦值過程，而“怎么做”即運算過程。一、矩陣回到之前介紹令窗口。開始學(xué)習(xí)如何賦值例 1 矩陣的直接賦值，生成矩

4、陣， = = %清除工作間中的變量，且命令必須小寫%清除顯示的內(nèi)容clearclc%在里的；表示換行X=2 1;1 -3Y=0 4;2 3X*YX.*Y運行結(jié)果：X =%矩陣乘法% 點乘，即矩陣對應(yīng)位置數(shù)相乘21其他略1-3計算 + = 例 2 = %將 X 的值賦予 A，；表示當(dāng)行命令不顯示%因為有 AX=B，所以 X=1B 可以用逆矩陣來計算% inv 表示對矩陣求逆A=X;B=13;-18;X=inv(A)*B*（或者）X=AB運行結(jié)果：X=37%這個是號，表示右邊的數(shù)除以左邊的數(shù)；%每一次變量被賦值的時候，它之前的內(nèi)容都會被替換掉例 3 增量賦值法輸入矩陣這里要教大家一種常用的冒號表達

5、式：%e1 為起點，e2 為步長，e3 為終點e1:e2:e3%自動生成從 1 到 10 等差數(shù)a=1:1:10運行結(jié)果：a =12345678910例 4 一個簡單的應(yīng)用，利用上面所學(xué)的冒號表達式生成下面矩陣12243615482051025612307891035404550： A=1:1:10; B=A;A*2;A*5對于冒號的用法還有C=B(3,3)C=B(:,2)D=B(3,:)大家還可以試下C=B(3)C=B(1:3,3 5)D=Bd= rot90(B,2)dragon=B(:)例 5 特殊矩陣的創(chuàng)建A=eye(3)B=zeros(3)C=ones(3,6)% B 矩陣中第三行第三

6、列的數(shù)：15%B 矩陣中第二列%B 矩陣中第三行%第三行第一列 5%從第一行到第三行的第三列和第五列%轉(zhuǎn)置%天翻地覆的感覺 。，逆時針旋轉(zhuǎn) 2*90 度%你猜猜.*C%從 1 到 10 生成 10 個等差數(shù)，與冒號的區(qū)別是不需要知道步長b=linspace(1,10,10)運行結(jié)果：b =12345678910%從 100 到 102 生成 11 個等比數(shù)d=logspace(0,2,11)d=1.00001.58492.51193.98116.309610.000015.848925.118939.810763.0957100.0000例 6 矩陣的結(jié)構(gòu)，練習(xí)將上述的 A 、B、 C 矩陣組

7、合起來%方括號中，逗號的意思是列UOU=A,B;C運行結(jié)果: UOU =100111whos010111001111000111000111000111%變量的詳細情況%返回矩陣的維度%可以試試a,b=size(UOU)a,b,c=size(UOU)clearwhos例 7 一個比較難的應(yīng)用，同學(xué)根據(jù)上面知識來生成YOY =100000：010000001000000111000111000111001000010000100000 A=zeros(3); B=ones(3); C=eye(3); D=rot90(C,1); YOY=C,A,D;A,B,A%將 C 矩陣逆時針旋轉(zhuǎn) 90 度二、

8、繪圖現(xiàn)在，利用剛才所學(xué)到令來嘗試做作圖x=0: 5:360*pi/180;plot(x,sin(x),x,cos(x),x,0);%plot 二維作圖中最基本運行結(jié)果如右圖所示：令例 8SIN(X)與 COS(X)作圖這里開始需要用到 M 文件，在講M 文件前下路徑問題，簡單點說就是，你需要將會用到的東西都放在一個文件夾里面，這樣才能順利的調(diào)用。這個路徑的修改可以在命令窗口里完成。點擊命令窗口左上角的圖標或快捷鍵“ctrl+N”新建一個clear; clf程序，編寫如下命令：NUMERIC=xlsread(sinx.xlsx) ;%從工作區(qū)域獲得名為 sinx.xlsx 中的數(shù)據(jù)，并賦值給變量

9、%第一列賦值給 x1%第二列賦值給 x2%用函數(shù)命令生成 x3x1=NUMERIC(:,1);x2=NUMERIC(:,2);x3=cos(x1);%當(dāng)然，以上步驟你也可以通過矩陣賦值替換。直接使用 x1=0: 5:360*pi/180;% or項的作用 o 是改變線形狀，r 是改變顏色plot(x1,x2,or,x1,x3,-g);%定義顯示坐標(xmin xmax ymin ymax)%圖的名字axis(0 6 -1 1)title(正余弦);xlabel(x);ylabel(sinx);作圖就學(xué)到這里，講義后面的附錄里有其他的作圖程序，有可以去看下。邏輯與循環(huán)有時需要判斷值之間的大小關(guān)系

10、，有時候需要對一種運算重復(fù)運行很多次，這些情況下，邏輯運算和循環(huán)語句就會被用到。邏輯運算有大小判斷(=；; i=(7=8)運行結(jié)果：i =0例 10 一個單循環(huán)，假設(shè) y=1/3+1/5+1/(2n-1)，求 y3 時的最大 n 值以及對應(yīng)的 y 值主程序：lizi_10.m建立clc,clear; y=0; n=1;while y3程序，輸入以下命令：%當(dāng) y=tol)k(i+1)=(1-alpha)*k(i)alpha)/(1+n)*(2+rho);i=i+1;endi;k_ss plot(1:i,k)運行結(jié)果是：k_ss =0.2052三、函數(shù)的調(diào)用1、背景：它是迭代法（切）：迭代法（N

11、ewtons method）又稱為-遜方法（Newton-Raphson method），在 17 世紀一種在實數(shù)域和復(fù)數(shù)域上近似求解方程的方法。該方法的思維：第一步，設(shè) r 是 f(x) = 0 的根，選取 x0 作為 r 的初始近似值，過點(x0，f(x0)做曲線 y = f(x)的切線 L，第二步，L 的方程為 y = f(x0)+f(x0)(x-x0），求出 L 與 x 軸交點的橫坐標 x1 = x0-f(x0)/f(x0)，稱 x1 為 r 的一次近似值。第三步，過點（x1，f(x1)做曲線 y = f(x)的切線，并求該切線與 x 軸交點的橫坐標 x2 = x1-f(x1)/f(x

12、1)，稱 x2 為 r 的二次近似值。第四步，重復(fù)以上過程，得 r 的近似值序列，其中，x(n+1)=x(n）f(x(n)/f(x(n)稱為 r 的 n+1次近似值。例 12 運用法求方程 = 的根在命令窗口輸入： syms x; f=(x3-2*x-1) x,k=Newtondd(f,-0.9,10-8)運行結(jié)果是：x =-1y =6%定義 x 為符號變量%定義你所求的方程%運用新建的函數(shù)求解2、初值問題的數(shù)值解法經(jīng)典-法（RK4）背景：-法（Runge-Kutta）是用于模擬常微分方程的解的重要的一類隱式或顯式迭代法。這些技術(shù)由數(shù)學(xué)家和馬丁海姆于 1900 年左右發(fā)明。初值問題表述為： y

13、/ t = f (t，y(t) ，已經(jīng)初始狀態(tài)y(0) = 0，用數(shù)值方法求出某時間段 t=0 tn內(nèi)在步長間隔上的各個時刻狀態(tài)變量 y（t）的數(shù)值解，即求解過程順著節(jié)點排列的次序一步一步地向前推進。則該問題的 RK4 由下面的方程給出：n+1 = + h/6(1 + 22 + 23 + 4)其中：1 = (, )2 = ( + 2 , + 2 1)3 = ( + 2 , + 2 2)4 = ( + , + 3)這樣，下一個值(n+1)由現(xiàn)在的值()加上時間間隔(h)和一個估算的斜率的乘積決定。例 13 用四階-法求下面常微分方程的數(shù)值解= + ( + )() = 在一個 M 文件中編寫微分方

14、程 function dy=q(t,y,m) dy=1+log(t+1);之后，在命令窗口輸入clear all( )%注意單引號，在 RK4 文件中時間矩陣被定義為列t=0:0.1:1 ; t,k=rk4(q,t,1)運行結(jié)果：k = 1.00001.10481.21881.34111.47111.60821.75201.90212.05802.21952.3863例 14 常微分方程組的數(shù)值解 = + + = ( ) + () = () = 在一個 M 文件中編寫微分方程組function dy=q_2(t,y,m) dy=y(1)+y(2)+2;y(1)-y(2)+2;之后，在命令窗口輸

15、入clear allt=0:0.1:1 ;%初值定義y0=2 1;t,k=rk4(q_2,t,y0)部分運行結(jié)果：t =00.10001.0000k =2.00002.54171.00001.311013.55426.2831Ramsey 模型的數(shù)值模擬回顧一下講義中的 Ramsey 模型，具體微分方程組為：= ( ) （ ） = ( + )參數(shù)如下： = . nl = . = . = . = 根據(jù)微分方程組求出穩(wěn)態(tài)解 k_ss、c_ss，并基于這些條件就有右上的相位圖。_ 那么問題來了，在給定初始資本 k(0) 的條件下，如何選擇最優(yōu)的初始消費 c(0)。打靶法（shooting metho

16、d）這里需要解釋一種叫做打靶法的方法，它是一種求解微分方程組初值問題的數(shù)值算法。該方法的思路：第一步，依據(jù)條件約束確定初始值可取值的區(qū)間c_l c_h。第二步，選擇取值區(qū)間里的中點值作為初始值c_half (即 1/2*(c_l+c_h)來進行推演。若結(jié)果偏大則初始點選擇過高，也即是說中點值以上的點都偏高，應(yīng)該舍去。在這種情況下，取值的可選區(qū)間重新改寫為c_l c_half，即原中點值成為新的取值上限 c_h。反之，若結(jié)果偏小，則舍去下半部分區(qū)間，取值的可選區(qū)間重新寫為c_half c_h，即原中點成為新的取值下限。第三步，基于新的取值區(qū)間再次選擇中點值作為初始值，并重復(fù)上面的步驟，直至選出滿

17、足精度的最優(yōu)初始點。運用這個算法，求解上述問題，具體的流程圖，寫在了下一頁。本題的函數(shù)文件是：ramsey_rk4_practice.m主程序文件是：shootingramsey_rk4_practice.m（Tip：在執(zhí)行程序前要記得先把工作路徑改到 mat 文件夾。詳情參考繪圖那一頁）Ramsey 的函數(shù)文件如下：function dx=ramsey_rk4_practice(t,x,flag)global A alpha1 nl rhx1=x(1); %消費 x2=x(2); %資本heta;dx1=(x1/theta)*(alpha1*x2(alpha1-1)-delta-rho);d

18、x2=x2alpha1-x1-(nl+delta)*x2; dx=dx1;dx2;流程圖：主程序：clearglobal alpha1 nl rh tic;heta;%計時插件% % *% parameters% * alpha1=0.3;nl=0.07;rho=0.02; delta=0.02; theta=3;%參數(shù)設(shè)置% *% steady-se values% * chi_s=(delta+rho)/alpha1)(-1/(1-alpha1);c_s=(delta+rho)/alpha1)(-alpha1/(1-alpha1)-(nl+delta)*(delta+rho)/alpha1

19、)(-1/(1-alpha1);%穩(wěn)態(tài)值disp(c_s,chi_s);% *% initial conditions% *%輸入資本初值chi_0=input(chi_0 = );% *% resolving IVPs% *% 時間的步長% 精度%消費額下限%消費額上限%選取消費中間值嘗試hh=0.1;tol=1e-4; c_l=0;c_h=chi_0alpha1;c_0=(c_l+c_h)/2;%給定初始條件X0=c_0,chi_0;%試驗計數(shù)器%間隔，第一次進入iter=1;crit=1;循環(huán)的條件%最大實驗次數(shù)%觸邊條件，如果 k(t) c(t)沒有達到穩(wěn)態(tài)，就繼續(xù)循環(huán)%第一次進入內(nèi)循

20、環(huán)所需條件maxit=200;while (abs(crit)=tol) dd=1;i=1; t=1;disp(X0); X=X0;while min(dd)=0%顯示初值情況%當(dāng)沒有出現(xiàn)遞減的情況下，繼續(xù)循環(huán)tn,xn=rk4(ramsey_rk4_practise,t t+hh,X0);dd=diff(xn);%利用 RK4 來計算下期的 c k%路徑X0=xn(2,:);X=X;xn(2,:);t=t+hh; i=i+1;end disp(dd); disp(X(i,1);crit=X(i,1)-c_s;disp(crit); disp(running time%顯示內(nèi)循環(huán)結(jié)束后的消費值

21、inutes: , num2str(toc/60), , after , num2str(iter), iterations);%如果消費出現(xiàn)下降的情況i(1)maxit) breakend%消費額重置，再次設(shè)定新的初始消費來嘗試%試驗次數(shù)計數(shù)enddisp(Convergence after , num2str(iter), iterations); disp(running time: , num2str(toc);%將消費值的路徑賦值給 x1%將資本值的路徑賦值給 x2% 返回 x1 的維度% 計算時間跨度x1=X(:,1);x2=X(:,2);a,b=size(x1); n=hh*a-

22、hh; t=0:hh:n;disp(c_s,chi_s);save callocation.mat X% *% results% *figure;plot(t,x1); grid on;title(optimal consumption path); xlabel(t);ylabel(c);figure;plot(t,x2); grid on;title(optimal capital path); xlabel(t);ylabel(chi);%開啟網(wǎng)格本節(jié)課的內(nèi)容就講到這里，以下是一些關(guān)于作圖，微分計算和邏輯判斷的例題，因為上間有限所以就沒能一一提到。附錄1 隱函數(shù)作圖：常常會遇見一些隱函數(shù)

23、作圖，這些函數(shù)大多是多值函數(shù)，常規(guī)的作圖法顯得很繁瑣，接著將介紹一點簡單的方法。對于f=(x,y)=0 的二維隱函數(shù)的圖像，可以直接用ezplot 函數(shù)制圖。繪出x2 + 2 = 10的圖像建立程序，輸入以下命令：%定義函數(shù) f=x2+y2-10%畫出 x2+y2-10=0 的圖像f=(x,y)x.2+y.2-10;ezplot(f) axis equal;運行結(jié)果是：2.球體作圖程序 clear allf=(x,y,z)x.2+y.2+z.2-10;x,y,z=meshgrid(linspace(-4,4,25); val=f(x,y,z);p,v=isosurface(x,y,z,val,

24、0);patch(fa,p,vertiview(3); grid on; axis equal運行結(jié)果：,v,facevertexcdata,jet(size(v,1),facecolor,w,edgecolor,flat);3.微分方程組的動態(tài)圖像程序題目為：解微分方程組1 = 2 32 = 1 33 = 0.51 1 21(0) = 0，2(0) = 1，3(0) = 1函數(shù)文件 rigid： function dy=rigid(t,y) dy(1)=y(2)*y(3);dy(2)=-y(1)*y(3);dy(3)=-0.51*y(1)*y(2);dy=dy(:);%對 dy 矩陣賦值M

25、文件：clear ally0=0 1 1;t,y=ode45(rigid,0 12,y0); plot(t,y(:,1),-,t,y(:,2),*,t,y(:,3),+)計算微分方程的特解 y=dsolve(D2y+4*Dy+29*y=0,y(0)=0,Dy(0)=15,x)%括號中第一部分是方程，中間是初值設(shè)定，最后是標明對 x 求導(dǎo)運行結(jié)果：y =(3*sin(5*x)/exp(2*x)y= /( + 定積分)在一個 M 文件中寫下：function f =fx(x)f=x.*sin(x)./(1+cos(x).*cos(x);%定義一個叫做 fx 的函數(shù)%需要注意的是這里是點乘、點除%保

26、存該文件，名字必須與函數(shù)名字一致：fx之后，在命令窗口輸入： y=quad(fx,0,pi)或者 y=quad(fx,0,pi)運行結(jié)果是：y =2.4674%quad 是積分令，后兩位是積分下限和上限%有時候，人們喜歡用這種形式的假如 a=4，b=9 判斷兩者的大小建立新a=4;b=9;i1=(ab)，輸入以下命令：% 定義ab的返還值% “”表示括號中的逆命題i3=(ab)在這里因為文件中并未定義函數(shù)function 所以文件名字可以隨便取，這樣不取到特殊意義的就行。其實做 M 文件。運行結(jié)果是：i1 =i3 =在中有含 function =f 之類的叫做函數(shù)文件，而其他叫%邏輯判斷返還值

27、中 1 表示“是”% 0 表示“否”10Ramsey 模型的另外一個程序（區(qū)別在于沒有使用 RK4 來求 t+1 期）建立程序，輸入如下命令：%清屏clear all; clf;A=2; beta=0.96; sigma=2; alpha=0.3;%設(shè)置參數(shù)k_ss=(beta*A*alpha)(1/(1-alpha);%建立 2 個邊際值c_ss=A*k_ssalpha-k_ss;%精度設(shè)置tol=0.001;% 獲得穩(wěn)態(tài)點左邊的穩(wěn)定路徑k0=0.01;cl=0; ch=A*k0alpha-k0; cm=(cl+ch)/2;%資本初始值%消費值上下限cl ch;cl 表示下限% ch 表示下

28、限%選擇一個中間值來試試%將初始資本賦予第一期資本變量 k1，將 cm 賦予第一期消費%時間變量設(shè)置k(1)=k0; c(1)=cm;t=1;%循環(huán)開始while (norm(k(t) c(t)-k_ss c_ss)c_ss) | (c(t)tol)| (norm(k(t) c(t)-k_ss c_ss)c_ss) & (norm(k(t) c(t)-k_ss c_ss)tol)%如果突破C 線則減小初始消費的上限ch=cm;elseif (c(t)tol) & (norm(k(t) c(t)-k_ss c_ss)tol)%如果突破K 線則增加初始消費的下限cl=cm;%else%break;%否則跳出該循環(huán)%里層循環(huán)的END%外層循環(huán)的ENDendendsaddle_path=k;c;% 獲得穩(wěn)態(tài)右邊的穩(wěn)定路徑clear k0 cl ch cm k c t; k0=0.9;cl=A*k0alpha-k0; ch=A*k0alpha