matlab第四讲习题
第四讲上机练习（1）用plot(x)命令画直线。
x1=[1 2 3]；x2=[0 1 0]（2）绘制正弦曲线y=sin(x)和方波曲线。
（3）利用矩阵绘制图线，观察结果。
X=peaks；plot（X）（4）在一个图形窗口内，利用plot 绘制多条三角函数曲线，注意默认的结果特征。
*hold on（5）用不同线段类型、颜色和数据点形画出sinx 和cosx 曲线。
（6）用图形表示连续调制波形y = sin(t)sin(9t) ，数据点之间不用曲线连接。
（7）将图形窗口分割为四个子图，并绘制适当的图形。
*subplot（8）在图形窗口中增加标识。
基于x=0: 0.1 : 2*pi ;plot(x,sin(3*x))；hold on；plot(x,cos(2*x),'ro')；增加图标题、坐标轴标识、图例并在特定位置增加文字注释。
*text、legend、title、xlabel（9）利用帮助学习plotyy 的功能，并利用它在同一图形窗口中绘制sinx 和2cos(2x) 的图形。
（10）三维曲线绘图，运行以下命令，注意观察结果。
x = 0 : 0.1 : 20*pi;plot3 ( x , sin(x) , cos(x) )（11）以sin(t) 、cos(t) 和cos(2*t) 为坐标绘制三维曲线。
（12）使用meshgrid命令并观察结果。
x = [ 1 2 3 4 ];y = [ 5 6 7 ];[ xx , yy ] = meshgrid( x , y )（13）测试meshgrid的使用方法。
x=linspace(-3,3,49);y=linspace(-3,3,49);[xx,yy]=meshgrid(x,y) ; %产生49*49的栅格点坐标mesh(xx) %查看xx的网线图mesh(yy)zz=3*(1-xx).^2.*exp(-(xx.^2) - (yy+1).^2) ...- 10*(xx/5 - xx.^3 - yy.^5).*exp(-xx.^2-yy.^2) ...- 1/3*exp(-(xx+1).^2 - yy.^2); %产生peaks函数plot3(xx,yy,zz)（14）上题中再使用mesh 、surf 命令绘制结果图。
（15）用meshz和meshc查看peaks 函数的三维曲面图。
*[X，Y，Z] = peaks(100)（16）使用以下修饰命令对图形进行调整。
colormap(MAP)colorbar()shading faceted / flat / interphidden on/offview(az,el)waterfall / meshzcontour(Z,n) / contour3(Z,n)（17）特殊命令绘图。
4.3 MATLAB的特殊图形绘制条形图常用于对统计的数据进行作图，特别适用于少量且离散的数据。
绘制条形图的函数如表4.10所示。
语法：bar(x,y,width,'参数') %画条形图bar3(y,z,width,'参数') %画三维条形图说明：x是横坐标向量，省略时默认值是1:m，m为y的向量长度；y是纵坐标，可以是向量或矩阵，当是向量时每个元素对应一个竖条，当是m×n的矩阵时，将画出m组竖条每组包含n条；width是竖条的宽度，省略时默认宽度是0.8，如果宽度大于1，则条与条之间将重叠；'参数'有grouped(分组式)和stacked(累加式)，省略时默认为grouped。
bar3命令的格式也相同，y必须是单调增加或减小，省略时为1:m；'参数'除了grouped和stacked还有detached(分离式)。
【例4.19】用条形图表示某年一月份中3日～6日连续四天的温度数据，y矩阵的各列分别表示平均温度、最高温度和最低温度，如图4.23所示，用条形图和三维条形图分别表示。
x=3:6;y=[5.3000 13.0000 0.40005.1000 11.8000 -1.70003.7000 8.1000 0.60001.5000 7.7000 -4.5000]bar(x,y) %画条形图bar3(x,y) %画三维条形图程序分析：由上图看出条形图是按行分组的，每组为每天的平均温度、最高温度和最低温度。
1. 面积图面积图是在曲线与横轴之间填充颜色，用于绘制面积图的命令为“area”，只能用于二维绘图。
语法：area(y) %画面积图area(x,y)说明：y可以是向量或矩阵，如果y是向量则绘制的曲线和plot命令相同，只是曲线和横轴之间填充颜色，如果y是矩阵则每列向量的数据构成面积叠加起来；x是横坐标，当x省略时则横坐标为1:size(y,1)。
2. 实心图实心图是将数据的起点和终点连成多边形，并填充颜色，绘制实心图的命令为“fill”。
语法：fill(x,y,c) %画实心图说明：c为实心图的颜色，可以用'r'、'g'、'b'、'c'、'm'、'y'、'w'、'k'，或RGB三元组行向量表示，也可以省略。
【例4.19续】绘制面积图和实心图，并比较其区别，如图4.24所示。
area(x,y) %面积图fill(x,y,'r') %红色的实心图程序分析：由上图可知面积图是绘制曲线和横轴间的面积，y的各列叠加在一起的，而实心图是将起点和终点连接并填充颜色的多边形。
语法：hist(y,m) %统计每段的元素个数并画出直方图hist(y,x)说明：m是分段的个数，省略时则默认为10；x是向量，用于指定所分每个数据段的中间值；y可以是向量或矩阵，如果是矩阵则按列分段。
【例4.20】用直方图表示正态分布的随机数分布，如图4.25所示。
y=randn(10,2) %产生10*2的正态分布的随机数矩阵y =-1.1878 -1.1859-2.2023 -1.05590.9863 1.4725-0.5186 0.05570.3274 -1.21730.2341 -0.04120.0215 -1.1283-1.0039 -1.3493-0.9471 -0.2611-0.3744 0.9535x=-2:0.5:2;hist(y,x)程序分析：直方图显示的是y在x附近的元素的个数，如-2附近有一个。
产生的随机数不同则得出的直方图也不同。
饼图是用于显示向量中的各元素占向量元素总和的百分比，可以用pie和pie3命令分别绘制二维和三维饼图。
语法：pie(x,explode,’label’) %画二维饼图pie3(x,explode,’label’) %画三维饼图说明：x是向量；explode是与x同长度的向量，用来决定是否从饼图中分离对应的一部分块，非零元素表示该部分需要分离；’label’是用来标注饼图的字符串数组。
【例4.21】绘制四个季度支出额的饼图，如图4.26所示。
y=[200 100 250 400]; %四个季度支出额explode=[0 0 1 0];pie(y,explode,{'第一季度','第二季度','第三季度','第四季度'})MATLAB 提供了多个绘制离散数据的命令，有stem 、stem3、stairs和scatter 等。
【例4.22】使用几种绘制离散数据的命令来显示)sin(2x e y x -=的离散数据，如图4.27所示。
x=0:0.1:2*pi; y=sin(x).*exp(-2*x); subplot(3,1,1) stem(x,y,'filled') %画火柴杆图 subplot(3,1,2) stairs(x,y) %画阶梯图 subplot(3,1,3)scatter(x,y) %画点图程序分析：'filled'参数是来填充火柴杆图的点标记。
1. 对数坐标图形对数坐标图形有semilogx 、semilogy 和loglog 命令。
semilogx(x,y,'参数') %绘制x 为对数坐标的曲线 semilogy(x,y,'参数') %绘制y 为对数坐标的曲线loglog(x,y,'参数') %绘制x 、y 都为对数坐标的曲线说明：参数和plot 命令一样，只是坐标不同。
【例4.23】求传递函数为)1s 5.0(s 1)s (G +=的对数幅频特性曲线，如图4.28所示，横坐标为w 按对数坐标。
w=logspace(-2,3,20); %频率w 为0.01到1000 Aw=1./(w.*sqrt((0.5*w).^2+1)); %计算幅频 Lw=20*log10(Aw); %计算对数幅频 semilogx(w,Lw)title('对数幅频特性曲线')2. 极坐标图极坐标图由polar 命令来实现。
polar(theta,radius,'参数') %绘制极坐标图说明：theta为相角，radius为离原点的距离。
【例4.23续】用极坐标图表示上述传递函数的Nyquist曲线，如图4.29所示。
w=logspace(-2,3,20);Fw=-90-atan(0.5*w);polar(Fw,Aw)语法：contour(Z,n) %绘制Z矩阵的等高线contour(x,y,z,n) %绘制以x和y指定x、y坐标的等高线说明：n为等高线的条数，省略时为自动条数。
【例4.24】绘制peaks函数的等高线，如图4.30所示。
[x,y,z]=peaks;contour(x,y,z) %画二维等高线contour3(z,30) %画30条三维等高线1. compass命令compass绘制的是以原点为起点的一组复向量，因此又称为罗盘图。
语法：compass(u,v) %画罗盘图compass(Z)说明：u、v分别为复向量的实部和虚部；当只有一个参数Z时，则相当于compass(real(Z),imag(Z))。
2. feather命令feather绘制的是起点为(k,0)的复向量图，又称为羽毛图。
语法：feather(u,v) %画羽毛图feather (Z)【例4.25】用罗盘图和羽毛图绘制复向量，如图4.31所示。
theta=0:0.2:2*pi;z=sin(theta).*exp(j*theta);compass(z)feather(z)程序分析：羽毛图的绘制起点是(k,0)，k从1～n，n是Z向量的元素序号。
第四讲 MATLAB绘图
希腊字母、上标、下标、数学符号、字型：
\ alpha \ beta \ gamma \ pi \ tall
\ Delta
\ delta
\ Omega
a2 a^{2} a2 a _{2} \ inf ty \ times \ oplus \ otimes
t = -pi:pi/100:pi; y = sin(t); plot(t,y) axis([-pi pi -1 1]) xlabel('-\pi \leq {\itt} \leq \pi’, 'FontSize',16) ylabel('sin(t) ', 'FontSize',16) title('Graph of the sine function') text(1,-1/3,'{\itNote the odd symmetry.}')
plot(x1, y1, 选项1, x2, y2, 选项2, …, xn, yn, 选项n) plot (x, y, ‘color_linestyle_marker’) 例 ： plot (x, y, ‘y:square’)
color_linestyle_marker
Color strings are 'c', 'm', 'y', 'r', 'g', 'b', 'w', and 'k'. These correspond to cyan, magen, white, and black.
%加图形标题
xlabel('independent variable X');
《MATLAB及其工程应用》第四讲
Plotting
How great if I can create a beautiful graph with MATLAB!
It is really not so difficult as you imaged.
What’s the first ?
• Right! You need a piece of paper. figure(number)
‘fprintf’ will only display the real part of a complex number!!! example: a = 3+4i; fprintf(‘fprintf: a = %f’, a);
Summary
There are 3 ways to display a data: 1) Leave the semicolon off the end of the statement 2) Use ‘disp’ function 3) Use ‘fprintf’ function It is better to use ‘disp’ to display a complex number.
disp
disp( string)
str = [ ‘the value of pi = ‘ num2str(pi) ]; disp(str);
fprintf
fprintf( format, data)
fprintf(‘the value of pi is %f \n’, pi);
Programming pitfalls
Data Formats
• • • • • • • • • • • • Short Long Short e Short g Long e Long g Bank Hex Rat Compact Loose +
Matlab学习教程 第四章(4)上机练习
第4章图形处理功能1 内容简介基本内容主要包括：（1）二维图形（2）三维图形（3）图形处理的基本技术2 达到的目标（1）掌握二维图形的绘制。
（2）掌握三维图形的绘制。
（3）掌握图形处理的基本技术3 具体内容３．1 二维图形３．1．1 基本绘图命令（1）当plot函数仅有一个输入变量例4-1y=[5 2 3 8 5]; %y 行矩阵plot(y) %一条线例4-2y=[5 2 3 8 5;2 4 3 1 5;1 1 1 1 1]; %y 矩阵plot(y) %5条线,等于矩阵的列数（2）当plot函数有两个输人变量例4-3x=0:0.01*pi:pi;y=sin(x).*cos(x);plot(x,y)例4-3x=0:0.01*pi:pi;y=[sin(x);cos(x); sin(x).*cos(x)];plot(x,y)例4-4x1=0:0.01*pi:pi;x2=pi:0.01*pi:2*pi;x=[x1' x2'];y=[sin(x1') cos(x2')];plot(x,y)例4-5x1=1:5;x2=6:10;y1=x1;y2=2*x2;plot([x1;x2],[y1;y2])%plot([x1' x2'],[y1' y2'])（3）当plot函数有三个输入变量时MATLAB语言中提供的对曲线的线型、颜色以及标识的控制符如表4．l所示。
例4-6 绘制带有显示属性设置的二维图形。
x=0.5*pi: 0.1*pi:2*pi;y=sin(x);z=cos(x);plot (x, y, '--ko', x, z, '-. r*')３．1．2 特殊的二维图形函数（1）特殊坐标系的二维图形函数（ａ）对数坐标例4-7 绘制X坐标为对数坐标的二维图形。
x=0.5*pi: 0.1*pi:2*pi;y=sin(x);semilogx (x, y, '-ro')（ｂ）极坐标例4-8绘制极坐标下的二维图形。
matlab第四章课后作业解答
matlab第四章课后作业解答第四章习题解答1、求下列多项式的所有根，并进行验算。
（3）267235865x x x x-+-（4）4)32(3-+x 解：>> p=zeros(1,24); >> p(1)=5;p(17)=-6;p(18)=8;p(22)=-5; >> root=roots(p)root =0.97680.9388 + 0.2682i0.9388 - 0.2682i0.8554 + 0.5363i0.8554 - 0.5363i0.6615 + 0.8064i0.6615 - 0.8064i0.3516 + 0.9878i0.3516 - 0.9878i-0.0345 + 1.0150i-0.0345 - 1.0150i-0.4609 + 0.9458i-0.4609 - 0.9458i-0.1150 + 0.8340i-0.1150 - 0.8340i-0.7821 + 0.7376i-0.7821 - 0.7376i-0.9859 + 0.4106i-0.9859 - 0.4106i-1.0416-0.7927>> polyval(p,root)ans =1.0e-012 *-0.07120.0459 - 0.0081i0.0459 + 0.0081i-0.0419 + 0.0444i-0.0419 - 0.0444i0.0509 + 0.0929i0.0509 - 0.0929i-0.2059 + 0.0009i-0.2059 - 0.0009i-0.0340 + 0.0145i-0.0340 - 0.0145i0.1342 + 0.0910i0.1342 - 0.0910i0.0025 + 0.0027i0.0025 - 0.0027i-0.0077 + 0.4643i-0.0077 - 0.4643i-0.3548 - 0.1466i-0.3548 + 0.1466i-0.0251-0.0073(4) >> p1=[2 3];>> p=conv(conv(p1,p1),p1)-[0 0 0 4]; >> root=roots(p)root =-1.8969 + 0.6874i-1.8969 - 0.6874i-0.7063>> polyval(p,root)ans =1.0e-014 *-0.7105 - 0.6217i-0.7105 + 0.6217i6、求解下列方程组在区域1,0<<βα内的解-=+=.sin 2.0cos 7.0,cos 2.0sin 7.0βαββαα 解：以初值)5.0,5.0(),(00=βα进行求解>> fun=inline('[0.7*sin(x(1))+0.2*cos(x(2))-x(1),0.7*cos(x(1))-0.2*sin(x(2))-x(2)]');>> [x,f,h]=fsolve(fun,[0.5 0.5])Optimization terminated: first-order optimality is less than options.TolFun.x =0.5265 0.5079f =1.0e-007 *-0.1680 -0.2712h =1因而，该方程组的近似根为5079.0,5265.0==βα。
matlab第4讲
2013-7-9
Matlab Language
14
6、算术运算 （续）
2013-7-9
Matlab Language
15
6、算术运算 （续）
【例5-2】点幂“.^”举 例 >>a=1:6
a= 1 2 3 4 5 6
>>a=a.^2
a= 1 4 9 16 25 36
>>b=b.^2
b= 1 4
>>b=reshape(a,2,3)
2013-7-9
Matlab Language
22
【例7-1】求向量的最大值 >>x=[-43,72,9,16,23,47]; >>y=max(x) %求向量x中的最大值 y= 72 >>[y,l]=max(x) %求向量x中的最大值及其该元素的位置 y= 72 l= 2
2013-7-9
Matlab Language
7
5、多维数组 （续）
三维数组元素的寻址：可以(行、列、页)来确定。 以维数为 3×4×2 的三维数组为例，其寻址方式如 下图所示：
数组 A 是三维数组，其中 A(：,：,1)代表第一页的二 维数组，A(：,：,2)代表第二页的二维数组。
Matlab Language
8
2013-7-9
5、多维数组 （续）
标量关系进行比较，并给出结果，形成一个维数与原来相同
的0、1矩阵。 3、当一个标量与一个矩阵比较时，该标量与矩阵的各元素进行
比较，结果形成一个与矩阵维数相等的0、1矩阵。
2013-7-9
Matlab Language
17
7、关系运算 （续） 【例】建立5阶方阵A，判断其元素能否被3整除。
第四讲matlab插值、拟合和回归分析
第四讲matlab插值、拟合和回归分析第四讲插值、拟合与回归分析在⽣产实践和科学研究中，常常有这样的问题：由实验或测量得到变量间的⼀批离散样本点，要求得到变量之间的函数关系或得到样本点之外的数据。
解决此类问题的⽅法⼀般有插值、拟合和回归分析等。
设有⼀组实验数据0011(,),(,),(,)n n x y x y x y ，当原始数据精度较⾼，要求确定⼀个简单函数()y x ?=（⼀般为多项式或分段多项式）通过各数据点，即(),0,,i i y x i n ?== ，称为插值问题。
另⼀类是拟合问题，当我们已经有了函数关系的类型，⽽其中参数未知或原始数据有误差时，我们确定的初等函数()y x ?=并不要求经过数据点，⽽是要求在某种距离度量下总体误差达到最⼩，即(),0,,i i i y x i n ?ε=+= ，且20ni i ε=∑达到最⼩值。
对同⼀组实验数据，可以作出各种类型的拟合曲线，但拟合效果有好有坏，需要进⾏有效性的统计检验，这类问题称为回归分析。
⼀、插值(interpolation)常⽤的插值⽅法有分段线性插值、分段⽴⽅插值、样条插值等。
1、⼀元插值yi=interp1(x,y,xi,method)对给定数据点(x,y)，按method 指定的⽅法求出插值函数在点（或数组）xi 处的函数值yi 。
其中method 是字符串表达式，可以是以下形式：'nearest' ——最邻近点插值'linear' ——分段线性插值（也是缺省形式）'spline' ——分段三次样条插值'cubic' 分段⽴⽅插值例：在⼀天24⼩时内，从零点开始每间隔2⼩时测得环境温度数据分别为（℃）：12，9，9，10，18，24，28，27，25，20，18，15，13⽤不同的插值⽅法估计中午1点（即13点）的温度，并绘出温度变化曲线。
>> x=0:2:24;>> y=[12 9 9 10 18 24 28 27 25 20 18 15 13];>>y_linear=interp1(x,y,13),y_nearest=interp1(x,y,13,'nearest')>>y_cubic=interp1(x,y,13,'cubic'),y_spline=interp1(x,y,13,'spline')>> y1=interp1(x,y,xx); y2=interp1(x,y,xx,'nearest');>> y3=interp1(x,y,xx,'cubic');y4=interp1(x,y,xx,'spline');>> subplot(2,2,1),plot(x,y,'or',xx,y1)>> subplot(2,2,2),plot(x,y,'or',xx,y2)>> subplot(2,2,3),plot(x,y,'or',xx,y3)>> subplot(2,2,4),plot(x,y,'or',xx,y4)2、⼆元插值zi=interp2(X,Y,Z,xi,yi,method)已知数据点（X，Y，Z），求插值函数在（xi,yi）处的函数值zi,插值⽅法method同interp1。
Matlab求解微分方程(组)及偏微分方程(组)
第四讲 Matlab 求解微分方程（组）理论介绍：Matlab 求解微分方程（组）命令 求解实例：Matlab 求解微分方程（组）实例实际应用问题通过数学建模所归纳得到的方程，绝大多数都是微分方程，真正能得到代数方程的机会很少.另一方面，能够求解的微分方程也是十分有限的，特别是高阶方程和偏微分方程（组）.这就要求我们必须研究微分方程（组）的解法：解析解法和数值解法. 一．相关函数、命令及简介1.在Matlab 中，用大写字母D 表示导数，Dy 表示y 关于自变量的一阶导数，D2y 表示y 关于自变量的二阶导数，依此类推.函数dsolve 用来解决常微分方程（组）的求解问题，调用格式为：X=dsolve(‘eqn1’,’eqn2’,…)函数dsolve 用来解符号常微分方程、方程组，如果没有初始条件，则求出通解，如果有初始条件，则求出特解.注意，系统缺省的自变量为t2.函数dsolve 求解的是常微分方程的精确解法，也称为常微分方程的符号解.但是，有大量的常微分方程虽然从理论上讲，其解是存在的，但我们却无法求出其解析解，此时，我们需要寻求方程的数值解，在求常微分方程数值解方面，MATLAB 具有丰富的函数，我们将其统称为solver ，其一般格式为：[T,Y]=solver(odefun,tspan,y0)说明：(1)solver 为命令ode45、ode23、ode113、ode15s 、ode23s 、ode23t 、ode23tb 、ode15i 之一.(2)odefun 是显示微分方程'(,)y f t y =在积分区间tspan 0[,]f t t =上从0t 到f t 用初始条件0y 求解.(3)如果要获得微分方程问题在其他指定时间点012,,,,f t t t t 上的解，则令tspan 012[,,,]f t t t t =(要单调的).(4)因为没有一种算法可以有效的解决所有的ODE 问题，为此，Matlab 提供了多种求解器solver ，对于不同的ODE 问题，采用不同的solver.表1 Matlab中文本文件读写函数说明：ode23、ode45是极其常用的用来求解非刚性的标准形式的一阶微分方程（组）的初值问题的解的Matlab常用程序，其中：ode23采用龙格-库塔2阶算法，用3阶公式作误差估计来调节步长，具有低等的精度.ode45则采用龙格-库塔4阶算法，用5阶公式作误差估计来调节步长，具有中等的精度.3．在matlab命令窗口、程序或函数中创建局部函数时，可用联函数inline，inline函数形式相当于编写M函数文件，但不需编写M-文件就可以描述出某种数学关系.调用inline函数，只能由一个matlab表达式组成，并且只能返回一个变量，不允许[u,v]这种向量形式.因而，任何要求逻辑运算或乘法运算以求得最终结果的场合，都不能应用inline函数，inline函数的一般形式为：FunctionName=inline(‘函数容’, ‘所有自变量列表’)例如：（求解F(x)=x^2*cos(a*x)-b ,a,b是标量；x是向量）在命令窗口输入:Fofx=inline(‘x .^2*cos(a*x)-b’ , ‘x’,’a’,’b’); g= Fofx([pi/3 pi/3.5],4,1) 系统输出为：g=-1.5483 -1.7259注意：由于使用联对象函数inline 不需要另外建立m 文件，所有使用比较方便，另外在使用ode45函数的时候，定义函数往往需要编辑一个m 文件来单独定义，这样不便于管理文件，这里可以使用inline 来定义函数. 二．实例介绍1.几个可以直接用Matlab 求微分方程精确解的实例 例1 求解微分方程2'2x y xy xe -+=程序：syms x y; y=dsolve(‘Dy+2*x*y=x*exp(-x^2)’,’x ’)例 2 求微分方程'0x xy y e +-=在初始条件(1)2y e =下的特解并画出解函数的图形.程序：syms x y; y=dsolve(‘x*Dy+y-exp(1)=0’,’y(1)=2*exp(1)’,’x ’);ezplot(y)例 3 求解微分方程组530tdx x y e dtdy x y dt⎧++=⎪⎪⎨⎪--=⎪⎩在初始条件00|1,|0t t x y ====下的特解并画出解函数的图形.程序：syms x y t[x,y]=dsolve('Dx+5*x+y=exp(t)','Dy-x-3*y=0','x(0)=1','y(0)=0','t') simple(x); simple(y)ezplot(x,y,[0,1.3]);axis auto2.用ode23、ode45等求解非刚性标准形式的一阶微分方程（组）的初值问题的数值解（近似解）例 4 求解微分方程初值问题2222(0)1dy y x xdx y ⎧=-++⎪⎨⎪=⎩的数值解，求解围为区间[0,0.5].程序：fun=inline('-2*y+2*x^2+2*x','x','y');[x,y]=ode23(fun,[0,0.5],1); plot(x,y,'o-')例 5 求解微分方程22'2(1)0,(0)1,(0)0d y dyy y y y dt dtμ--+===的解，并画出解的图形.分析：这是一个二阶非线性方程，我们可以通过变换，将二阶方程化为一阶方程组求解.令12,,7dyx y x dtμ===，则 121221212,(0)17(1),(0)0dx x x dtdx x x x x dt⎧==⎪⎪⎨⎪=--=⎪⎩ 编写M-文件vdp.m function fy=vdp(t,x)fy=[x(2);7*(1-x(1)^2)*x(2)-x(1)]; end在Matlab 命令窗口编写程序 y0=[1;0][t,x]=ode45(vdp,[0,40],y0);或[t,x]=ode45('vdp',[0,40],y0); y=x(:,1);dy=x(:,2); plot(t,y,t,dy)练习与思考：M-文件vdp.m 改写成inline 函数程序？ 3.用Euler 折线法求解Euler 折线法求解的基本思想是将微分方程初值问题00(,)()dyf x y dxy x y ⎧=⎪⎨⎪=⎩ 化成一个代数(差分)方程，主要步骤是用差商()()y x h y x h +-替代微商dydx，于是00()()(,())()k k k k y x h y x f x y x h y y x +-⎧=⎪⎨⎪=⎩记1,(),k k k k x x h y y x +=+=从而1(),k k y y x h +=+于是0011(),,0,1,2,,1(,).k k k k k k y y x x x h k n y y hf x y ++=⎧⎪=+=-⎨⎪=+⎩例 6 用Euler 折线法求解微分方程初值问题22(0)1dyx y dxy y ⎧=+⎪⎨⎪=⎩的数值解（步长h 取0.4），求解围为区间[0,2].分析：本问题的差分方程为00110,1,0.4,0,1,2,,1(,).k k k k k k x y h x x h k n y y hf x y ++===⎧⎪=+=-⎨⎪=+⎩程序：>> clear >> f=sym('y+2*x/y^2'); >> a=0; >> b=2; >> h=0.4; >> n=(b-a)/h+1; >> x=0; >> y=1;>> szj=[x,y];%数值解 >> for i=1:n-1y=y+h*subs(f,{'x','y'},{x,y});%subs ，替换函数 x=x+h; szj=[szj;x,y]; end >>szj>> plot(szj(:,1),szj(:,2))说明：替换函数subs 例如：输入subs(a+b,a,4) 意思就是把a 用4替换掉，返回 4+b ，也可以替换多个变量，例如：subs(cos(a)+sin(b),{a,b},[sym('alpha'),2])分别用字符alpha 替换a 和2替换b ，返回 cos(alpha)+sin(2)特别说明：本问题可进一步利用四阶Runge-Kutta 法求解，Euler 折线法实际上就是一阶Runge-Kutta 法，Runge-Kutta 法的迭代公式为001112341213243(),,(22),6(,),0,1,2,,1(,),22(,),22(,).k k k k k k k k k k k k y y x x x h h y y L L L L L f x y k n h h L f x y L h h L f x y L L f x h y hL ++=⎧⎪=+⎪⎪=++++⎪⎪=⎪=-⎨⎪=++⎪⎪⎪=++⎪⎪=++⎩相应的Matlab 程序为：>> clear >> f=sym('y+2*x/y^2'); >> a=0; >> b=2; >> h=0.4; >> n=(b-a)/h+1; >> x=0; >> y=1;>> szj=[x,y];%数值解 >> for i=1:n-1l1=subs(f, {'x','y'},{x,y});替换函数 l2=subs(f, {'x','y'},{x+h/2,y+l1*h/2}); l3=subs(f, {'x','y'},{x+h/2,y+l2*h/2}); l4=subs(f, {'x','y'},{x+h,y+l3*h}); y=y+h*(l1+2*l2+2*l3+l4)/6; x=x+h; szj=[szj;x,y]; end>>szj>> plot(szj(:,1),szj(:,2))练习与思考：(1)ode45求解问题并比较差异. (2)利用Matlab 求微分方程(4)(3)''20y y y -+=的解.(3)求解微分方程''2',2(1)0,030,(0)1,(0)0y y y y x y y --+=≤≤==的特解. (4)利用Matlab 求微分方程初值问题2''''00(1)2,|1,|3x x x y xy y y ==+===的解. 提醒：尽可能多的考虑解法 三．微分方程转换为一阶显式微分方程组Matlab 微分方程解算器只能求解标准形式的一阶显式微分方程（组）问题，因此在使用ODE 解算器之前，我们需要做的第一步，也是最重要的一步就是借助状态变量将微分方程（组）化成Matlab 可接受的标准形式.当然，如果ODEs 由一个或多个高阶微分方程给出，则我们应先将它变换成一阶显式常微分方程组.下面我们以两个高阶微分方程组构成的ODEs 为例介绍如何将它变换成一个一阶显式微分方程组.Step 1 将微分方程的最高阶变量移到等式左边，其它移到右边，并按阶次从低到高排列.形式为：()'''(1)'''(1)()'''(1)'''(1)(,,,,,,,,,,)(,,,,,,,,,,)m m n n m n x f t x x x x y y y y y g t x x x x y y y y ----⎧=⎨=⎩Step 2 为每一阶微分式选择状态变量，最高阶除外'''(1)123'''(1)123,,,,,,,,,m m n m m m m n x x x x x x x x x y x y x y x y--++++========注意：ODEs 中所有是因变量的最高阶次之和就是需要的状态变量的个数，最高阶的微分式不需要给它状态变量.Step 3 根据选用的状态变量，写出所有状态变量的一阶微分表达式''''122334123''12123,,,,(,,,,,),,(,,,,,)m m n m m m nm n x x x x x x x f t x x x x xx xg t x x x x +++++======练习与思考：(1)求解微分方程组**'''3312*'''3312()()22x x x y x r r y y y x y r r μμμμμμ⎧+-=+--⎪⎪⎨⎪=+--⎪⎩其中2r =1r =*1,μμ=-1/82.45,μ=(0) 1.2,x =(0)0,y ='(0)0,x ='(0) 1.049355751y =-(2)求解隐式微分方程组''''''''''''2235x y x y x y x y xy y ⎧+=⎨++-=⎩ 提示：使用符号计算函数solve 求'''',x y ，然后利用求解微分方程的方法 四．偏微分方程解法Matlab 提供了两种方法解决PDE 问题，一是使用pdepe 函数，它可以求解一般的PDEs,具有较大的通用性，但只支持命令形式调用；二是使用PDE 工具箱，可以求解特殊PDE 问题，PDEtoll 有较大的局限性，比如只能求解二阶PDE 问题，并且不能解决片微分方程组，但是它提供了GUI 界面，从复杂的编程中解脱出来，同时还可以通过File —>Save As 直接生成M 代码.1.一般偏微分方程（组）的求解(1)Matlab 提供的pdepe 函数，可以直接求解一般偏微分方程（组），它的调用格式为：sol=pdepe(m,pdefun,pdeic,pdebc,x,t)pdefun 是PDE 的问题描述函数，它必须换成标准形式：(,,)[(,,,)](,,,)m m u u u uc x t x x f x t u s x t u x t x x x-∂∂∂∂∂=+∂∂∂∂∂ 这样，PDE 就可以编写入口函数：[c,f,s]=pdefun(x,t,u,du)，m,x,t 对应于式中相关参数，du 是u 的一阶导数，由给定的输入变量可表示出c,f,s 这三个函数.pdebc 是PDE 的边界条件描述函数，它必须化为形式：(,,)(,,).*(,,,)0up x t u q x t u f x t u x∂==∂ 于是边值条件可以编写函数描述为：[pa,qa,pb,qb]=pdebc(x,t,u,du),其中a 表示下边界，b 表示上边界.pdeic 是PDE 的初值条件，必须化为形式：00(,)u x t u =，故可以使用函数描述为：u0=pdeic(x)sol 是一个三维数组，sol(:,:,i)表示i u 的解，换句话说，k u 对应x(i)和t(j)时的解为sol(i,j,k)，通过sol ，我们可以使用pdeval 函数直接计算某个点的函数值.(2)实例说明 求解偏微分2111222221220.024()0.17()u u F u u t xu u F u u tx ⎧∂∂=--⎪⎪∂∂⎨∂∂⎪=+-⎪∂∂⎩ 其中， 5.7311.46()x x F x e e -=-且满足初始条件12(,0)1,(,0)0u x u x ==及边界条件1(0,)0,u t x ∂=∂221(0,)0,(1,)1,(1,)0uu t u t t x∂===∂ 解：(1)对照给出的偏微分方程和pdepe 函数求解的标准形式，原方程改写为111221220.024()1.*()10.17u u F u u x u F u u u t x x ∂⎡⎤⎢⎥--⎡⎤⎡⎤⎡⎤∂∂∂=+⎢⎥⎢⎥⎢⎥⎢⎥-∂∂∂⎣⎦⎣⎦⎣⎦⎢⎥⎢⎥∂⎣⎦可见1121220.024()10,,,()10.17u F u u x m c f s F u u u x ∂⎡⎤⎢⎥--⎡⎤⎡⎤∂====⎢⎥⎢⎥⎢⎥-∂⎣⎦⎣⎦⎢⎥⎢⎥∂⎣⎦ %目标PDE 函数function [c,f,s]=pdefun(x,t,u,du) c=[1;1];f=[0.024*du(1);0.17*du(2)]; temp=u(1)-u(2);s=[-1;1].*(exp(5.73*temp)-exp(-11.46*temp)) end(2)边界条件改写为：下边界2010.*00f u ⎡⎤⎡⎤⎡⎤+=⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦上边界1110.*000u f -⎡⎤⎡⎤⎡⎤+=⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦%边界条件函数function [pa,qa,pb,qb]=pdebc(xa,ua,xb,ub,t) pa=[0;ua(2)]; qa=[1;0]; pb=[ub(1)-1;0]; qb=[0;1]; end(3)初值条件改写为：1210u u ⎡⎤⎡⎤=⎢⎥⎢⎥⎣⎦⎣⎦%初值条件函数 function u0=pdeic(x) u0=[1;0]; end(4)编写主调函数 clc x=0:0.05:1; t=0:0.05:2; m=0;sol=pdepe(m,pdefun,pdeic,pdebc,x,t); subplot(2,1,1) surf(x,t,sol(:,:,1)) subplot(2,1,2) surf(x,t,sol(:,:,2))练习与思考： This example illustrates the straightforward formulation, computation, and plotting of the solution of a single PDE.2()u u t x xπ∂∂∂=∂∂∂ This equation holds on an interval 01x ≤≤ for times 0t ≥. The PDE satisfies the initial condition (,0)sin u x x π= and boundary conditions(0,)0;(1,)0t uu t e t xπ-∂=+=∂ 2.PDEtool 求解偏微分方程(1)PDEtool （GUI ）求解偏微分方程的一般步骤在Matlab 命令窗口输入pdetool ，回车，PDE 工具箱的图形用户界面(GUI)系统就启动了.从定义一个偏微分方程问题到完成解偏微分方程的定解，整个过程大致可以分为六个阶段Step 1 “Draw 模式”绘制平面有界区域Ω，通过公式把Matlab 系统提供的实体模型：矩形、圆、椭圆和多边形，组合起来，生成需要的平面区域.Step 2 “Boundary 模式”定义边界，声明不同边界段的边界条件.Step 3 “PDE 模式”定义偏微分方程，确定方程类型和方程系数c,a,f,d ，根据具体情况，还可以在不同子区域声明不同系数.Step 4 “Mesh 模式”网格化区域Ω，可以控制自动生成网格的参数，对生成的网格进行多次细化，使网格分割更细更合理.Step 5 “Solve 模式”解偏微分方程，对于椭圆型方程可以激活并控制非线性自适应解题器来处理非线性方程；对于抛物线型方程和双曲型方程，设置初始边界条件后可以求出给定时刻t 的解；对于特征值问题，可以求出给定区间上的特征值.求解完成后，可以返回到Step 4，对网格进一步细化，进行再次求解.Step 6 “View 模式”计算结果的可视化，可以通过设置系统提供的对话框，显示所求的解的表面图、网格图、等高线图和箭头梯形图.对于抛物线型和双曲线型问题的解还可以进行动画演示.(2)实例说明用法求解一个正方形区域上的特征值问题：12|0u u u u λ∂Ω⎧-∆-=⎪⎨⎪=⎩ 正方形区域为：11,1 1.x x -≤≤-≤≤(1)使用PDE 工具箱打开GUI 求解方程(2)进入Draw 模式，绘制一个矩形，然后双击矩形，在弹出的对话框中设置Left=-1,Bottom=-1,Width=2,Height=2,确认并关闭对话框(3)进入Boundary 模式，边界条件采用Dirichlet 条件的默认值(4)进入PDE 模式，单击工具栏PDE 按钮，在弹出的对话框中方程类型选择Eigenmodes,参数设置c=1,a=-1/2,d=1,确认后关闭对话框(5)单击工具栏的 按钮，对正方形区域进行初始网格剖分，然后再对网格进一步细化剖分一次(6)点开solve菜单，单击Parameters选项，在弹出的对话框中设置特征值区域为[-20,20](7)单击Plot菜单的Parameters项，在弹出的对话框中选中Color、Height(3-D plot)和show mesh项，然后单击Done确认(8)单击工具栏的“=”按钮，开始求解。
第4讲 MATLAB作图
在区间[0，10*pi]画出参数曲线x=sin(t),y=cos(t), z=t. 解 t=0:pi/50:10*pi; plot3(sin(t),cos(t),t) Matlab liti8 rotate3d %旋转
2、多条曲线 、 PLOT3(x,y,z)
其中x，y，z是都是m*n矩阵，其对应的每一列表示一条曲线. 例 画多条曲线观察函数Z=(X+Y).^2. 解 x=-3:0.1:3;y=1:0.1:5; [X,Y]=meshgrid(x,y); Z=(X+Y).^2; plot3(X,Y,Z) Matlab liti9
例 在区间[0,2*pi]画sin(x)的图形，并加注图例“自变量 X”、“函数Y”、“示意图”, 并加格栅. 解 x=linspace(0,2*pi,30); y=sin(x); plot(x,y) xlabel('自变量X') ylabel(' ylabel('函数Y') Y') title('示意图') grid on Matlab liti2
(3)meshz(X,Y,Z) 在网格周围画一个curtain图(如,参考平面) 例 绘peaks的网格图
解 输入命令: [X,Y]=meshgrid(-3:.125:3); Z=peaks(X,Y); Meshz(X,Y,Z) Matlab liti36
返回
图 形 处 理
在图形上加格栅、 在图形上加格栅、图例和标注 定制坐标 图形保持 分割窗口 缩放图形 改变视角 动 画
例 在[-1，2]上画 y = e
解
2x
+ sin(3x 2 ) 的 图形
Matlab liti43
第四讲 MATLAB与EXCEL数据交互
2 读取数据xlsread函数
xlsread函数语法 1．[ data,textdate]= xlsread(filename) 输入参数： Filename：目标文件地址（若文件在matlab当前 的工作目录中，Filename为’文件名’，如果文件 不在matlab当前的工作目录中，filename为’文件 路径\文件名’） 输出参数： Data: 数值数据 Textdate: 文字数据
2．data= xlsread(filename, sheet, raቤተ መጻሕፍቲ ባይዱge) 输入参数：
Filename：目标文件地址（若文件在matlab当前的 工作目录中，Filename为’文件名’，如果文件 不在matlab当前的工作目录中，filename为’文 件路径\文件名’）
Sheet：数据表名称，例如excel默认表名称sheet1。 Range：数据所在位置，例如A1，B13等 输出参数：
4.3 Excel2007加载与使用宏 加载方法： 点击excel的office按钮点击excel选项在加载项中 点击 转到见下图
浏览（matlab的安装路径）toolbox文件夹 exlink文件夹 excllink.xla文件（打开）
使用方法： excel2007加载项下可以发现exlink相关的按钮， 具体使用方法与exlink在excel2003中的使用方 法一样。
目录
中，filename为’文件路径\文件名’）； M: 写入excel中的数据，M存储数据的变量名称； Sheet: 写入excel中的sheet名称（ 可选，若空默认sheet1）； Range：写入excel中的单元格区域（可选，若空默认’A1’）； 输出参数： status: 写入状态 “1”表示写入成功“0”表示写入失败 message: 若失败，则显现失败信息
matlab第四讲课件共27页文档
格式化输出—fprintf
cows=5; fprintf('There are %f cows in the pasture',cows) There are 5.000000 cows in the pasture
格式化输出—fprintf
格式化输出—fprintf
Matlab在执行完函数fprintf后不会自动重起 一行。前述的命令行执行完后，在命令窗 口中命令提示符紧跟在函数输出字符串的 后面，并没有另起一行。若再次执行其它 命令，则两次的输出结果会在同一行中显 示出来。
这种函数格式将下列字符串写入文件my_output_file
中。
Some example output is 3141.29
并且在命令窗口返回写入数据的字节数： ans=32
注意
在使用fprintf时，初学者常犯的错误是 பைடு நூலகம்记在占位符后输入域类型标识，如f，这 样函数将不会正常工作，而且还不会给出 错误提示。
其中，函数fopen的第一个输入参数是要打开的文件
名。第二个输入参数是字符串‘wt’,表示要对文件进行
写的操作。如果能够正确打开这个输出文件，并且已经给
该文件分配了文件标识符，就可以把这个文件标识符作为
函数fprintf的第一个输入参数按照指定的格式把数据写
入到文件中：
fprintf(file_id,’some example output is %4.2f \n’,pi*1000)
cows=1:5; fprintf('There are %f cows in the pasture',cows)
格式化输出—fprintf
如果需要分行显示，则在字符串后使用\n进 行换行。