宏观课件matlab课堂程序例题

1、本节课的主要目的在于让了解一些的基本知识。因为大分部同学都没有怎么接触过这个，所以前半节我会讲一些简单的操作，大家跟着做下，基本上就能明白。而后半节课内容会难一点，但我也会尽力让大家明白的有趣之处。方面：的安装与运行。没有的同学可以下课来找我，然后我默认大家都会安装这个，并且能顺利打开它。打开之后，首先看见的界面叫做“命令窗口”，光标闪的地方就是你编辑命令的地方，在这个窗口每一句输入的指令都会立即执行。（尝试下 clc）而点击界面左上角第一个图标可以新建一个程序，界面像一张空白页，如下图所示，编入这里令不会立即执行，并且可以让你修改，所以这里才是程序编辑的地方。编程完成之后，你需要保存该文件，

2、操作是 FileSave as选择保存路径取个好名字哦保存。这样一来你就生成了一个后缀名为.m 的文件，这类文件就被叫做M 文件。之后，如何重新打开它们，例如打开 XX.m。操作是 FileOpen选择路径XX.m(XX:刚才你输入的那个文件名)。最后，在子程序的运行操作中点击图标或者快捷键”F5”便可运行程序，结果将会出现在之前令窗口中。的函数和命令有很多，大家想深入了解，可以有以下方法：如果你知道指令，但不清楚用法，可以使用自带的 help 命令 如：help global如果你不知道指令，那可以进行模糊搜索，在命令窗口输入相关词按 Tab 键就能出现所有相关令；3.当然你也可以“加相关指

3、令” 如：“clf”。学= =# 如果你连指令是什么都不是很清楚的话，就先把这下吧。（讲义中，表示在命令窗口输入令，而%表示的是注释，*后面的表示有难度可以跳过）的特点高效方便的矩阵数组运算灵活简单的绘图功能库函数丰富，避免了繁杂的子程序编程任务*众多面向领域应用的工具箱主要介绍的是前面的 3 个特点，而最后的工具箱主要面对较高水平的应用，所以本课是不会涉及到的。编辑程序：其实编程就是人机交流，计算机有自己能识别的语言，所以你得用这些计算机语言来告诉它是什么，怎么做。 这里，“是什么”即赋值过程，而“怎么做”即运算过程。一、矩阵回到之前介绍令窗口。开始学习如何赋值例 1 矩阵的直接赋值，生成矩

4、阵， = = %清除工作间中的变量，且命令必须小写%清除显示的内容clearclc%在里的；表示换行X=2 1;1 -3Y=0 4;2 3X*YX.*Y运行结果：X =%矩阵乘法% 点乘，即矩阵对应位置数相乘21其他略1-3计算 + = 例 2 = %将 X 的值赋予 A，；表示当行命令不显示%因为有 AX=B，所以 X=1B 可以用逆矩阵来计算% inv 表示对矩阵求逆A=X;B=13;-18;X=inv(A)*B*（或者）X=AB运行结果：X=37%这个是号，表示右边的数除以左边的数；%每一次变量被赋值的时候，它之前的内容都会被替换掉例 3 增量赋值法输入矩阵这里要教大家一种常用的冒号表达

5、式：%e1 为起点，e2 为步长，e3 为终点e1:e2:e3%自动生成从 1 到 10 等差数a=1:1:10运行结果：a =12345678910例 4 一个简单的应用，利用上面所学的冒号表达式生成下面矩阵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 特殊矩阵的创建A=eye(3)B=zeros(3)C=ones(3,6)% B 矩阵中第三行第三

6、列的数：15%B 矩阵中第二列%B 矩阵中第三行%第三行第一列 5%从第一行到第三行的第三列和第五列%转置%天翻地覆的感觉 。，逆时针旋转 2*90 度%你猜猜.*C%从 1 到 10 生成 10 个等差数，与冒号的区别是不需要知道步长b=linspace(1,10,10)运行结果：b =12345678910%从 100 到 102 生成 11 个等比数d=logspace(0,2,11)d=1.00001.58492.51193.98116.309610.000015.848925.118939.810763.0957100.0000例 6 矩阵的结构，练习将上述的 A 、B、 C 矩阵组

7、合起来%方括号中，逗号的意思是列UOU=A,B;C运行结果: UOU =100111whos010111001111000111000111000111%变量的详细情况%返回矩阵的维度%可以试试a,b=size(UOU)a,b,c=size(UOU)clearwhos例 7 一个比较难的应用，同学根据上面知识来生成YOY =100000：010000001000000111000111000111001000010000100000 A=zeros(3); B=ones(3); C=eye(3); D=rot90(C,1); YOY=C,A,D;A,B,A%将 C 矩阵逆时针旋转 90 度二、

8、绘图现在，利用刚才所学到令来尝试做作图x=0: 5:360*pi/180;plot(x,sin(x),x,cos(x),x,0);%plot 二维作图中最基本运行结果如右图所示：令例 8SIN(X)与 COS(X)作图这里开始需要用到 M 文件，在讲M 文件前下路径问题，简单点说就是，你需要将会用到的东西都放在一个文件夹里面，这样才能顺利的调用。这个路径的修改可以在命令窗口里完成。点击命令窗口左上角的图标或快捷键“ctrl+N”新建一个clear; clf程序，编写如下命令：NUMERIC=xlsread(sinx.xlsx) ;%从工作区域获得名为 sinx.xlsx 中的数据，并赋值给变量

9、%第一列赋值给 x1%第二列赋值给 x2%用函数命令生成 x3x1=NUMERIC(:,1);x2=NUMERIC(:,2);x3=cos(x1);%当然，以上步骤你也可以通过矩阵赋值替换。直接使用 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);作图就学到这里，讲义后面的附录里有其他的作图程序，有可以去看下。逻辑与循环有时需要判断值之间的大小关系

10、，有时候需要对一种运算重复运行很多次，这些情况下，逻辑运算和循环语句就会被用到。逻辑运算有大小判断(=；; i=(7=8)运行结果：i =0例 10 一个单循环，假设 y=1/3+1/5+1/(2n-1)，求 y3 时的最大 n 值以及对应的 y 值主程序：lizi_10.m建立clc,clear; y=0; n=1;while y3程序，输入以下命令：%当 y=tol)k(i+1)=(1-alpha)*k(i)alpha)/(1+n)*(2+rho);i=i+1;endi;k_ss plot(1:i,k)运行结果是：k_ss =0.2052三、函数的调用1、背景：它是迭代法（切）：迭代法（N

11、ewtons method）又称为-逊方法（Newton-Raphson method），在 17 世纪一种在实数域和复数域上近似求解方程的方法。该方法的思维：第一步，设 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 的二次近似值。第四步，重复以上过程，得 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)运行结果是：x =-1y =6%定义 x 为符号变量%定义你所求的方程%运用新建的函数求解2、初值问题的数值解法经典-法（RK4）背景：-法（Runge-Kutta）是用于模拟常微分方程的解的重要的一类隐式或显式迭代法。这些技术由数学家和马丁海姆于 1900 年左右发明。初值问题表述为： y

13、/ t = f (t，y(t) ，已经初始状态y(0) = 0，用数值方法求出某时间段 t=0 tn内在步长间隔上的各个时刻状态变量 y（t）的数值解，即求解过程顺着节点排列的次序一步一步地向前推进。则该问题的 RK4 由下面的方程给出：n+1 = + h/6(1 + 22 + 23 + 4)其中：1 = (, )2 = ( + 2 , + 2 1)3 = ( + 2 , + 2 2)4 = ( + , + 3)这样，下一个值(n+1)由现在的值()加上时间间隔(h)和一个估算的斜率的乘积决定。例 13 用四阶-法求下面常微分方程的数值解= + ( + )() = 在一个 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)运行结果：k = 1.00001.10481.21881.34111.47111.60821.75201.90212.05802.21952.3863例 14 常微分方程组的数值解 = + + = ( ) + () = () = 在一个 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)部分运行结果：t =00.10001.0000k =2.00002.54171.00001.311013.55426.2831Ramsey 模型的数值模拟回顾一下讲义中的 Ramsey 模型，具体微分方程组为：= ( ) （ ） = ( + )参数如下： = . nl = . = . = . = 根据微分方程组求出稳态解 k_ss、c_ss，并基于这些条件就有右上的相位图。_ 那么问题来了，在给定初始资本 k(0) 的条件下，如何选择最优的初始消费 c(0)。打靶法（shooting metho

16、d）这里需要解释一种叫做打靶法的方法，它是一种求解微分方程组初值问题的数值算法。该方法的思路：第一步，依据条件约束确定初始值可取值的区间c_l c_h。第二步，选择取值区间里的中点值作为初始值c_half (即 1/2*(c_l+c_h)来进行推演。若结果偏大则初始点选择过高，也即是说中点值以上的点都偏高，应该舍去。在这种情况下，取值的可选区间重新改写为c_l c_half，即原中点值成为新的取值上限 c_h。反之，若结果偏小，则舍去下半部分区间，取值的可选区间重新写为c_half c_h，即原中点成为新的取值下限。第三步，基于新的取值区间再次选择中点值作为初始值，并重复上面的步骤，直至选出满

17、足精度的最优初始点。运用这个算法，求解上述问题，具体的流程图，写在了下一页。本题的函数文件是：ramsey_rk4_practice.m主程序文件是：shootingramsey_rk4_practice.m（Tip：在执行程序前要记得先把工作路径改到 mat 文件夹。详情参考绘图那一页）Ramsey 的函数文件如下：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;%参数设置% *% 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);%稳态值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;%试验计数器%间隔，第一次进入iter=1;crit=1;循环的条件%最大实验次数%触边条件，如果 k(t) c(t)没有达到稳态，就继续循环%第一次进入内循

20、环所需条件maxit=200;while (abs(crit)=tol) dd=1;i=1; t=1;disp(X0); X=X0;while min(dd)=0%显示初值情况%当没有出现递减的情况下，继续循环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%显示内循环结束后的消费值

21、inutes: , num2str(toc/60), , after , num2str(iter), iterations);%如果消费出现下降的情况i(1)maxit) breakend%消费额重置，再次设定新的初始消费来尝试%试验次数计数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);%开启网格本节课的内容就讲到这里，以下是一些关于作图，微分计算和逻辑判断的例题，因为上间有限所以就没能一一提到。附录1 隐函数作图：常常会遇见一些隐函数

23、作图，这些函数大多是多值函数，常规的作图法显得很繁琐，接着将介绍一点简单的方法。对于f=(x,y)=0 的二维隐函数的图像，可以直接用ezplot 函数制图。绘出x2 + 2 = 10的图像建立程序，输入以下命令：%定义函数 f=x2+y2-10%画出 x2+y2-10=0 的图像f=(x,y)x.2+y.2-10;ezplot(f) axis equal;运行结果是：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运行结果：,v,facevertexcdata,jet(size(v,1),facecolor,w,edgecolor,flat);3.微分方程组的动态图像程序题目为：解微分方程组1 = 2 32 = 1 33 = 0.51 1 21(0) = 0，2(0) = 1，3(0) = 1函数文件 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)%括号中第一部分是方程，中间是初值设定，最后是标明对 x 求导运行结果：y =(3*sin(5*x)/exp(2*x)y= /( + 定积分)在一个 M 文件中写下：function f =fx(x)f=x.*sin(x)./(1+cos(x).*cos(x);%定义一个叫做 fx 的函数%需要注意的是这里是点乘、点除%保

26、存该文件，名字必须与函数名字一致：fx之后，在命令窗口输入： y=quad(fx,0,pi)或者 y=quad(fx,0,pi)运行结果是：y =2.4674%quad 是积分令，后两位是积分下限和上限%有时候，人们喜欢用这种形式的假如 a=4，b=9 判断两者的大小建立新a=4;b=9;i1=(ab)，输入以下命令：% 定义ab的返还值% “”表示括号中的逆命题i3=(ab)在这里因为文件中并未定义函数function 所以文件名字可以随便取，这样不取到特殊意义的就行。其实做 M 文件。运行结果是：i1 =i3 =在中有含 function =f 之类的叫做函数文件，而其他叫%逻辑判断返还值

27、中 1 表示“是”% 0 表示“否”10Ramsey 模型的另外一个程序（区别在于没有使用 RK4 来求 t+1 期）建立程序，输入如下命令：%清屏clear all; clf;A=2; beta=0.96; sigma=2; alpha=0.3;%设置参数k_ss=(beta*A*alpha)(1/(1-alpha);%建立 2 个边际值c_ss=A*k_ssalpha-k_ss;%精度设置tol=0.001;% 获得稳态点左边的稳定路径k0=0.01;cl=0; ch=A*k0alpha-k0; cm=(cl+ch)/2;%资本初始值%消费值上下限cl ch;cl 表示下限% ch 表示下

28、限%选择一个中间值来试试%将初始资本赋予第一期资本变量 k1，将 cm 赋予第一期消费%时间变量设置k(1)=k0; c(1)=cm;t=1;%循环开始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;%否则跳出该循环%里层循环的END%外层循环的ENDendendsaddle_path=k;c;% 获得稳态右边的稳定路径clear k0 cl ch cm k c t; k0=0.9;cl=A*k0alpha-k0; ch=A*k0alpha