matlab求matlab 解微分方程程my"+cy'+ky=F(t)，m,c,k均为常数，当F(t)以离散的表格形式给
点击联系发帖人
时间：2017-09-23 04:08
关于matlab解微分方程的问题，函数含有未知常数，怎么定义函数？_matlab吧_百度贴吧
&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&签到排名：今日本吧第个签到，本吧因你更精彩，明天继续来努力！
本吧签到人数：0成为超级会员，使用一键签到本月漏签0次！成为超级会员，赠送8张补签卡连续签到：天&&累计签到：天超级会员单次开通12个月以上，赠送连续签到卡3张
关注：172,543贴子：
关于matlab解微分方程的问题，函数含有未知常数，怎么定义函数？收藏
请问各位大神，如果我想把这个方程F(q) = dq/dt = q(1 - q)( pr4 - c2) F( p) = dp/dt = p(1 - p)(qr1 + qc3 - qr2 - c1)在matlab里面定义出来，然后再画出它的相位图，其中p,q是未知变量，其他的r1,r2,r4,c1,c2,c3是未知常数，最后求出的解中可以用其作为表达式，我该怎么写函数，写的函数好像都有问题，有大神帮忙吗？我写了两种，好像都通不过，大神帮忙看一下：一种：function dxdt=eq2(t,x)global a b c d e fdxdt=[x(1)*(1-x(1))*((a-b+e)*x(2)-d);x(2)*(1-x(2))*(c*x(1)-f)];end另一种：function xdot=eq3(t,x)%创建数组储存数据global a b c d e fxdot=zeros(2,1);xdot(1)=x(1)*(1-x(1))*((a-b+e)*x(2)-d);xdot(2)=x(2)*(1-x(2))*(c*x(1)-f);
网站建设10年沉淀，搭建专属品牌网站
楼主我也遇到这个问题了，微分方程里面有未知常数，这种方程怎么求解
登录百度帐号推荐应用查看: 30172|回复: 18|关注: 1
matlab dsolve求解微分方程组
<h1 style="color:# 麦片财富积分
新手, 积分 25, 距离下一级还需 25 积分
关注者： 5
最近有人问有关微分方程求解问题 特举一例供大家参考：
dsolve调用格式：dsolve（'equ1','equ2',..........'equN'）
另外要注意：在微分方程表达式输入中，以大写字母D来表示微分，D,D2,.......Dn分别表示一阶，二阶和n阶
2dx/dt+dy/dt-y=exp(-t)
dx/dt+x+y=0& && && && && && && && && &&&其中初始条件：x(0)=1.5,y(0)=0
首先求解微分方程的通解：
s=dsolve('2*Dx+Dy-y=exp(-t)','Dx+x+y=0');%求解微分方程组的通解
%微分方程组变量x的通解
-C1*exp((1+2^(1/2))*t)-C2*exp(-(2^(1/2)-1)*t)+1/2*C1*exp((1+2^(1/2))*t)*2^(1/2)-1/2*C2*exp(-(2^(1/2)-1)*t)*2^(1/2)-1/2*exp(-t)
%微分方程组变量y的通解
C1*exp((1+2^(1/2))*t)+C2*exp(-(2^(1/2)-1)*t)
然后根据初始条件，求解微分方程组的特解：
&& s=dsolve('2*Dx+Dy-y=exp(-t)','dx+x+y=0','x(0)=1.5','y(0)=0');%微分方程组在给定初始条件下的特解
Warning: Explicit solution could not be found.
& In dsolve at 333
&& s=dsolve('2*Dx+Dy-y=exp(-t)','Dx+x+y=0','x(0)=1.5','y(0)=0');%微分方程组在给定初始条件下的特解
-2^(1/2)*exp((1+2^(1/2))*t)+2^(1/2)*exp(-(2^(1/2)-1)*t)+exp((1+2^(1/2))*t)+exp(-(2^(1/2)-1)*t)-1/2*exp(-t)
2^(1/2)*exp((1+2^(1/2))*t)-2^(1/2)*exp(-(2^(1/2)-1)*t)
&& %或者使用下面命令直接获取x,y的特解
&& [x,y]=dsolve('2*Dx+Dy-y=exp(-t)','Dx+x+y=0','x(0)=1.5','y(0)=0')
-2^(1/2)*exp((1+2^(1/2))*t)+2^(1/2)*exp(-(2^(1/2)-1)*t)+exp((1+2^(1/2))*t)+exp(-(2^(1/2)-1)*t)-1/2*exp(-t)
2^(1/2)*exp((1+2^(1/2))*t)-2^(1/2)*exp(-(2^(1/2)-1)*t)
[ 本帖最后由 edifiers2008 于
11:12 编辑 ]
<h1 style="color:# 麦片财富积分
论坛优秀回答者
帖子最佳答案
关注者： 595
如果用simulink怎么求解非线性微分方程组呢?
<h1 style="color:# 麦片财富积分
<h1 style="color:# 麦片财富积分
请教您一个问题，为什么我的MATLAB RA2009中的dsolve函数不能调用呢？使用help dsolve命令说没有找到，用dsolve求解时就会报错说没有定义的函数。谢谢
<h1 style="color:# 麦片财富积分
我的、dsolve命令也不可用也是2009a
<h1 style="color:# 麦片财富积分
自己在网上找一个，拷进去就可以用
<h1 style="color:# 麦片财富积分
谢谢楼主！学习了，哈哈！
<h1 style="color:# 麦片财富积分
回复 6# damaomao 的帖子
2009a的dsolve命令不可用是怎么解决的？谢谢
<h1 style="color:# 麦片财富积分
新手，学习一下！
站长推荐 /2
快速搭建新能源汽车整车模型及其性能优化
MATLAB中文论坛是全球最大的 MATLAB & Simulink 中文社区。用户免费注册会员后，即可下载代码，讨论问题，请教资深用户及结识书籍作者。立即注册加入我们吧！
MATLAB官方社交平台
MATLAB中文论坛微社区4084人阅读
Matlab（10）
自动控制原理（3）
&[time,x]=solver(str,t,x0)&& 计算ODE或由字符串str 给定的ODE的值，部分解已在向量time中给出。在向量time
中给出部分解，包含的是时间值。还有部分解在矩阵x中给出，x的列向量每
个方程在这些值下的解。对于标量问题，方程的解将在向量x中给出。这些解
在时间区间t(1)到t(2)上计算得到。其初始值是x0 即x(t(1)).此方程组有str 指定的
M文件中函数表示出。这个函数需要两个参数：标量t 和向量x, 应该返回向量
x’(即x的导数)。因为对标量ODE来说，x和x’都是标量。在M文件中输入odefile
可得到更多信息。同时可以用命令numjac来计算Jacobi函数。
[t,x]=solver(str,t,x0,val)&& 此方程的求解过程同上，结构val包含用户给solver的命令。参见odeset和表1，可得到更多信息。
Ode45&&&&&此方法被推荐为首选方法。
Ode23&&&&&这是一个比ode45低阶的方法。
Ode113&&&&&用于更高阶或大的标量计算。
Ode23t&&&&&用于解决难度适中的问题。
Ode23s&&&&&用于解决难度较大的微分方程组。对于系统中存在常量矩阵的情况也有用。
Ode15s&&&&&与ode23相同，但要求的精度更高。
Ode23tb&&&&&用于解决难度较大的问题，对于系统中存在常量矩阵的情况也有用。
Set=odeset(set1,vak1,set2,val2,…)& 返回结构set,其中包含用于ODE求解方程的设置参数，&& 有关可用设置
的信息参见表1。
Odeget(set,’set1’)&&& 返回结构set中设置set1的值。
有许多设置对odeset控制的ODE解是有用的，参见表1。例如，如果要在求解过程中画出解的图形，可以
输入：inst=odeset(‘outputfcn’,’odeplot’);.也可使用命令odedemo。
表1&&& ODE求解方程的设置参数
RelTol&&&给出求解方程允许的相对误差
AbsTol&&给出求解方程允许的绝对误差
Refine&&给出与输入点数相乘的因子
OutputFcn&这是一个带有输入函数名的字符串，该字符串将在求解函数执行的每步被调用：odephas2(画出2D
的平面相位图)。Odephas3( 画出3D 的平面相位图)，odeplot( 画出解的图形)，odeprint( 显示中间结
OutputSel&是一个整型向量。指出哪些元素应该被传递给函数，特别是传递给OutputFcn
Stats&&如果参数Stats为on,则将统计并显示出计算过程中资源消耗情况
Jacobian… 如果编写ODE文件代码以便F（t,y,jocobian）返回dF/dy,则将jacobian设置为on
Jconstant…如果雅可比数df/dy是常量，则将此参数设置为on
Jpattern…& 如果编写ODE文件的编码以便函数F（[],[],jpattern）返回带有零的稀疏矩阵并输出非零元素dF/fy,则
需将Jpattern设置为on
Vectorized…如果编写ODE文件的编码以便函数F(t,[y1,y2……])返回[F(t,y1)& F(t,y2)…],则将此参数设置成on
Events…&& 如果ODE文件中带有参数‘events’，则将此参数设置成on
Mass…&& 如果编写ODE文件编码以实现函数F(t,[],‘mass’)返回M和M(t),应将此参数设置成on
MassConstant…如果矩阵M(t)是常量，则将此参数设置成on
MaxStep…&& 此参数是限定算法能使用的区间长度上限的标量
InitialStep…&& 给出初始步长的标量。如果给定的区间太大，算法就使用一个较小的步长
MaxOrder…&& 此参数只能被ode15s使用，它主要是指定ode15s的最高阶数,并且此参数应是从1到5的整数
BDF…&& 此参数只能被ode15s使用，如果倒推微分公式而不是用通常所使用的微分公式，则要将它设置为
NormControl…如果算法根据norm(e)&=max(Reltol*norm(y),Abstol)来步积分过程中的错误，则要将它设置成on
求解x’=-x^2 初值x(0)=1
创建函数xprim1,将此函数保存在M文件xprim1.m中：
function xprim=xprim1(t,x)
xprim=-x.^2;
然后调用MATLAB的ODE算法求解方程。然后画出解的图形：
[t,x]=ode45(‘xprim1’,[0,1],1);
plot(t,x,’-‘,t,x,’o’);
xlabel(‘time t0=0,tt=1’);
ylabel(‘x values x(0)=1’);
首先创建xprim2,将此函数保存在M文件xprim2.m中：
function xprim=xprim2(t,x)
xprim=x.^2;
然后调用MATLAB的ODE算法求解方程。然后画出解的图形：
[t,x]=ode45(‘xprim2’,[0,0.95],1);
注意：在MATLAB中计算出的点在微分绝对值大的区域内更密集些。
如果想在图像中直接观察，则在代码后输入以下命令
plot(t,x,’o‘,t,x,’-’);
xlabel(‘time t0=0,tt=0.95’);
ylabel(‘x values x(0)=1’);
例2求解下列二元一阶微分方程组
X1’=X1-0.1X1X2+0.01t
X2’=-X2+0.02X1X2+0.044
这个方程组用在人口动力学中。可以认为是单一化的捕食者---被捕食者模式。例如，狐狸和兔子。表示被捕食者， 表示捕食者。如果被捕食者有无限的食物，并且不会出现捕食者。于是有，这个式子是以指数形式增长的。大量的被捕食者将会使捕食者的数量增长；同样，越来越少的捕食者会使被捕食者的数量增长。而且，人口数量也会增长。
创建xprim3,将此函数保存在M文件xprim3.m 中：
function xprim=xprim3(t,x)
xprim=[x(1)-0.1*x(1)*x(2)+0.01*t; …
-x(2)+0.02*x(1)*x(2)+0.04*t];
然后调用一个ODE算法和画出解的图形：
[t,x]=ode45(‘xprim3’,[0,20],[30;20]);
plot(t,x);
xlabel(‘time t0=0,tt=20’);
ylabel(‘x values x1(0)=30,x2(0)=20’);
在MATLAB中,也可以根据x2函数绘制出x1图形,命令plot(x(:2),x(:1)) 可绘制出平面相位图.
例3 带参数的微分方程
X1’=a-(b+1)X1+X1^2X2
X2’=bX1+X1^2X2
方程由下面的M文件stiff1.m定义:
functionstiff=stiff1(t,x)
&& %变量不能放入参数表中
stiff=[0,0];&&& %stiff必须是一个冒号变量
stiff(1)=a-(b+1)*x1+x(1)^2*x(2);
stiff(2)=b*x(1)-x(1)^2*x(2);
下面的M文件给出一个比较困难的问题：
[t,X]=ode23(‘stiff1’，[0 10],[1 3]);
运行后得到的结果如下：
elapsed_time=72.1647
用专门解决复杂问题的解法ode23s,将得到较好的结果：
elapsed_time=1.0098
对于边界值问题，除了微分方程，还有边界处的值。在一维下这意味着至少有两个条件。
例4假设要研究一根杆的温度分布情况。这根杆一端的温度是T0,另一端的温度是T1。令y(x)表示这根杆的温度，函数f(x)表示加热源。从时间t=0开始，在相当长的时间内加热这根杆，直至达到平衡状态。这就是所谓的定常值或稳定状态。这个定常值可由下面的方程模型表示：
-Y’(X)=f(X),0&X&1
假设这根杆两端为：x=0和x=1。
假设在其两端又一根固定的柱子（或者可以看成是一个连接两个岛屿的桥），如图7所示。
令y(x)表示加载函数g(x)后弯曲的柱子。此问题需要有两个关于此柱子两端的边界条件。假设这根柱子非常牢
固的固定在墙上，即y在墙上的导数是0。可以得到下面的ODE,其中介绍了自然协调系统：
Y’’’(X)=g(x),0&X&1
Y(0)= Y’(0)= Y(1)= Y’(1)=0
解y(x)的导数可由有限的差分代替
如果用这些差分方程来代替ODE中的导数，就能得到一个所有未知的yi的方程组。其系数矩阵是一个有序区
间，此区间的宽度决定于这个微分方程的导数个数.
根据前面的温度模型的方程研究一下杆的温度分布，将所有的导数换成不同的差分并得到：
其中fj=f(xj)。为了简单起见，设M=6,即给定的y0和y6，而y1,y2,….y5 为未知变量。于是就有
注意，y0=T0和yM=T1必须移到方程组的右边。此时得到的矩阵是一个对角矩阵，其对角线上的元素为2，并且上一对角线和下一对角线上的元素为1。
下面解此问题的文件temperature.m。用户必须先给出分段数及f(x)(用点符号)，最后给出T0和T1。有关稀疏矩阵的更多信息参见其他资料。
%杆上的温度分布，用T0和T1分别表示两端温度
%这根杆放在x坐标的0和1 区间上，并被分成M个子区间，每个子区间的长度为1/M
%创建稀疏矩阵方程Ax=b并求解
%矩阵A对角阵，并以稀疏矩阵的形式存储
M=input(‘Give the number ofsubintervals (M):’);
Deltax=1/M;
xx=0:deltax:1;
funcStr= input(‘give f(x),the extra heatsource(e.g.,x.^3):’,’s’);
T0=input(‘Give y(0) (left): ’);
T1=input(‘Give y(1) (right): ’);
%构造对角阵和方程右边b
vectorOnes=ones(M-1,1);
A=spdiags([-vectorOnes,2*vectorOnes,-vectorOnes],[-1 0 1],M-1,M-1);
x=xx(2:end-1);
%x为区域内的值。
f=eval(funcStr);
%响应的f(x)的值。
b=deltax^2f;
b(1)=b(1)+T0; %对边界值x=0,x=1 进行特殊处理。
b(end)=b(end)+T1;
%解线性方程
%y在区间内：j=1:M-1.
y=[T0;y;T1];
%y在整个区间内：0&=x&=1.
%上面图形表示外部热源。
%下面图形表示杆上的热分布。
subplot(2,1,1);
plot(x,f);
title(‘External heatsource f(x).’,’FontSize’,14);
subplot(2,1,2);
plot(xx,y,’r’);
title(‘Tempearture distribution in a rod.’,’FontSize’,14);
%将区间分成等份，根据方程f(x)=在图中可以得到解。
例5 Matlab 中二阶常微分方程的数值解法
y''-3y'+2y=1;& y(0)=1;&y'(0)=0; 用数值解法求y(0.5)
odefun=@(t,x)[x(2);3*x(2)-2*x(1)+1];
[t,y]=ode45(odefun,[0:0.01:2],[1 0]);
y= dsolve('D2y-3*Dy+2*y=1','Dy(0)=0','y(0)=1');
y =exp(t) - exp(2*t)/2 + 1/2
&& feval(@(t)exp(t) - exp(2*t)/2 +1/2,0.5)
ans = 0.7896
例6求解范德普方程
首先是函数定义：
function xdot = myFcn (t, x)
xdot = [x(2) ; u(1-x(1)^2)*x(2)-x(1)]
其次是主程序：
% Solving options predetermined
options=odeset('OutputFcn',@odephas2,'reltol',1e-5);
% Solving vdp equation withdifferent ini conditions
for a = -3 : 0.1 : 3
b1=sqrt(25-a^2);
b2=-sqrt(25-a^2);
% Solve the problem
[t,y]=ode45(@myFcn, [0 20],[a b1], options);
[t,y]=ode45(@myFcn, [0 20],[a b2], options);
例7x1.^2-10x1+x2.^2+b=0
&&&& x1*x2.^2+x1-10x2+8=0
用MATLAB实现，编写非线性方程组的M文件fx1.m如下所示：
&& function y=fx1(x)
&&y(1)=x(1)*x(1)-10*x(1)+x(2)*x(2)+8;
&&y(2)=x(1)*x(2)*x(2)+x(1)-10*x(2)+8;
编写非线性方程组导数的M文件dfx1.m
& function y=dfx1(x)
&&&& y(1)=2*x(1)-10;
&&&& y(2)=2*x(2);
&&&& y(3)=x(2)*x(2)+1;
&&&& y(4)=2*x(1)*x(2)-10;
附一篇文章《用MATLAB求解非线性微分方程》
原始地址：http://blog.csdn.net/anhuixue/article/details/7560558
总结一下MATLAB中求解微分方程的思路和步骤。MATLAB中的帮助文件，内容实在是特别给力，简明的不能再简明了。可惜了是英文写的，即使是天天要面对英语的我，也颇为头疼。那么以下我将结合我对MATLAB中帮助文档的理解，提出我对求解最简单微分方程的方法的总结。以求解一个经典的非线性方程为例。
一求解一下这个方程，判断她是否稳定。要是稳定，那么她是否存在极限环;
二.一般的计算机求解方程的方法：
首先把该方程改写成一个规范的形式，一般使用状态空间表示法；
调用已有的算法进行求解；
最后对得出的结果进行处理，比如画图之类的。
三.输入待求解的方程。
MATLAB中关于输入输入待解方程的语句特别简单。需要先定义一个普通函数，函数名任意，姑且叫做myFcn，格式如下
function xdot = myFcn (t, x)&
需要注意的是，函数必须含有t, x两个参数，名称可以自己任意定。xdot是这个函数本身的返回值，只出现在这个函数内部，因此也可以任意定。
定义中的x被视为是一个列向量，x(i)表示列向量中的第i个分量。那么F函数的每一个分量即简单地用表达时给出即可。其中的自变量可以引用x(i)。以范德普方程为例：
xdot = [x(2) ; u(1-x(1)^2)*x(2)-x(1)]&
于是，这两句话便构成了待解函数。
四.调用MATLAB函数进行求解
通常人工求解微分方程需要知道初始值，计算机求解也不例外。另外，由于非线性方程一般只有数值解，故计算精度也可以调整。这些都是可以自己调整的参数。
调用MATLAB计算求解常微分方程的模式很简单，格式为：
[t, x] = ode45(@myFcn, [开始时间 结束时间], [初始值列向量],options)
注意到求解出来的结果是[t, x]即一堆数，所以需要我们进行后处理比如画图之类的。
ode45表示了MATLAB中的一种内置计算微分方程的算法，最为常用，出于娱乐目的，就先用这个。
@表示待求解的方程，如myFcn。
[开始时间 结束时间]表示求解的时间段，不必定义间隔。
[初始值列向量]就不用多说了。通常在求解之前定义一个变量x0列向量表示初始值，然后输入起来就方便多了。
options本身是一个变量。注意，她的名称不固定，但是他是一个以结构体为类型的变量。options很重要，稍后讨论。
五.进行后处理。在得到[t, x]后可以自己用plot神马的进行画图。
六.options的用法。
options在MATLAB帮助中是这样定义格式的：
&&&&&&options = odeset (“Name”, Value, “Name”, Value, …)
意思是，options这个变量要用到odeset这个函数来对他进行赋值。而odeset这个函数的参数必须是成对出现的：一个名称对应一个数值。
我们要用到的是这样的：
options= odeset(“OutputFcn”, @odephas2,“reltol”, 1e-5);
注意，odeset中的参数名称不可任意定，因为都是预定义好的。”OutputFcn”使用预定义的函数进行输出，而与定义的函数有好几种，使用@进行选择，我们选odephas2即画相平面图(phase portrait)。之后的那个参数表示相对误差的最大允许值，不说也明白。
有一个问题，用odephas2画出的图图不好看，因为上面打了好多圈圈。解决办法是找到该文件，修改其中所有的plot语句即可，把那个”-o”去掉就行了。
&&相关文章推荐
* 以上用户言论只代表其个人观点，不代表CSDN网站的观点或立场
访问：66045次
排名：千里之外
原创：28篇
转载：13篇
(3)(1)(1)(6)(9)(14)(7)
(window.slotbydup = window.slotbydup || []).push({
id: '4740881',
container: s,
size: '200,200',
display: 'inlay-fix'Matlab_学霸学习网
Matlab 软件目录1、Matlab 帮助的使用 ........................................................................ 3 2、Matlab 数据输入与类型 ................................................................ 4 3、Matlab 中的 M 文件及程序调试 ................................................ 12 4、Matlab 绘图命令 .......................................................................... 17 5、Matlab 在高等数学中的应用...................................................... 34 6、Matlab 在线性代数中的应用...................................................... 60 7、Matlab 数据处理 .......................................................................... 67 9、评价方法 ...................................................................................... 82 10、预测方法 .................................................................................... 97 11、蒙特卡洛方法 .......................................................................... 110 12、智能算法 .................................................................................. 122 13、分形 .......................................................................................... 129 14、Simulink 初步 .......................................................................... 134 15、Matlab 在概率统计中的应用.................................................. 147 参考文献 .......................................................................................... 1511 作为和 Mathematica、 Maple 并列的三大数学软件。 其强项就是其强大的矩阵计算以 及仿真能力。要知道 Matlab 的由来就是 Matrix + Laboratory = Matlab，所以这个软件在 国内也被称作《矩阵实验室》 。每次 MathWorks 发布 Matlab 的同时也会发布仿真工具 Simulink。在欧美很多大公司在将产品投入实际使用之前都会进行仿真试验，他们所主 要使用的仿真软件就是 Simulink。 Matlab 提供了自己的编译器： 全面兼容 C++以及 Fortran 两大语言。所以 Matlab 是工程师，科研工作者手上最好的语言，最好的工具和环境。 Matlab 已经成为广大科研人员的最值得信赖的助手和朋友！ 目前 MATLAB 产品族可以用来进行： － 数值分析 － 数值和符号计算 － 工程与科学绘图 － 控制系统的设计与方针 － 数字图像处理 － 数字信号处理 － 通讯系统设计与仿真 － 财务与金融工程... Simulink 是基于 MATLAB 的框图设计环境，可以用来对各种动态系统进行建模、 分析和仿真，它的建模范围广泛，可以针对任何能够用数学来描述的系统进行建模，例 如航空航天动力学系统、卫星控制制导系统、通讯系统、船舶及汽车等等，其中了包括 连续、离散，条件执行，事件驱动，单速率、多速率和混杂系统等等。 Simulink 提供 了利用鼠标拖放的方法建立系统框图模型的图形界面，而且 Simulink 还提供了丰富的 功能块以及不同的专业模块集合，利用 Simulink 几乎可以做到不书写一行代码完成整 个动态系统的建模工作。2 1、Matlab 帮助的使用1.1 helphelp?????????????%帮助总览 help elfun ???%关于基本函数的帮助信息 help exp ?????%指数函数 exp 的详细信息1.2 lookfor 指令当要查找具有某种功能但又不知道准确名字的指令时， help 的能力就不够了， lookfor 可以根据用户提供的完整或不完整的关键词，去搜索出一组与之相关的指令。 lookfor integral ??%查找有关积分的指令 lookfor fourier ??%查找能进行傅里叶变换的指令1.3 超文本格式的帮助文件在 Matlab 中，关于一个函数的帮助信息可以用 doc 命令以超文本的方式给出，如?doc ??? doc doc ??? doc eig ????%eig 求矩阵的特征值和特征向量1.4 pdf 帮助文件可从 MathWorks 网站上下载有关的 pdf 帮助文件。 网站地址：/3 2、Matlab 数据输入与类型2.1 Matlab 中的变量MATLAB 程序中的基本数据单元称为阵列(Array)，是一个分为行与列的数据集合。 变量被看作是只有一行一列的阵列。MATLAB 语言不需要对变量进行事先声明，也不 需要指定变量类型， 它会自动根据所赋予变量的值或对变量所进行的操作来确定变量的 类型。其命名规则为： （1）变量名的大小写是敏感的。 （2）变量的第一个字符必须为英文字母，而且不能超过 31 个字符。 （3）变量名可以包含下划线、数字，但不能为空格符、标点。 （4）命名变量时可以取一个容易记忆并且能表达出其含义的名称，如汇率，可以 定义为 exchange_rate。 对于变量作用域，默认情况是局部变量，使用 global 定义全局变量，而且全局变 量常用大写的英文字母表示。 MATLAB 预定义的变量如下表所示： ans eps pi inf NaN i 或 j nargin nargout realmax realmin flops 预设的计算结果的变量名 MATLAB 定义的正的极小值=2.2204e-16 内建的π 值 ∞值，无限大 无法定义一个数目 虚数单位 i=j=√-1 函数输入参数个数 函数输出参数个数 最大的正实数 2 最小的正实数 2 浮点运算次数41023-1022 注 1：在定义变量时要尽量与避免与这些名字相同，以免改变它们的值，如果已经 改变，可以通过 clear 复这些值。 注 2：数字的输入输出格式。所有数据均按 IEEE 浮点标准的长型格式存储。输入 格式沿用了 C 语言的风格和规则；输出格式使用 format 数据格式 命令控制，只影响 在屏幕上的显示结果，不影响内部的存储和运算。 变量名 来恢复它的初始值，也可以通过重新启动 MATLAB 恢2.2 向量及其运算1、向量的生成 ①命令窗口直接输入，使用[ ]，元素之间用空格、逗号或者分号隔开。 ②使用冒号表达式，基本形式为 x=x0:step:xn，其中 xn 为尾元素数值限，而不一定 是尾元素的值。当 step=1 时可省略步长。 ③生成线性等分向量，使用 linspace 函数。Y=linspace(x1,x2,n) ④生成对数等分向量，使用 logspace 函数。Y=logspace(x1,x2,n) 2、向量的基本运算5 数加（减） 、数乘、点积(dot 函数)、叉积(cross 函数)、混合积 dot(a,cross(b,c))2.3 矩阵及其运算1、简单矩阵的输入 （1）要直接输入矩阵时，矩阵一行中的元素用空格或逗号分隔；矩阵行与行之间 用分号“;”隔离，整个矩阵放在方括号“[ ]”里。 A=[1,2,3;4,5,6;7,8,9] ??? 说明：指令执行后，矩阵 A 被保存在 Matlab 的工作空间中，以备后用。如果用户 不用 clear 指令清除它，或对它进行重新赋值，那么该矩阵会一直保存在工作空间中， 直到本次指令窗关闭为止。 （2）矩阵的分行输入，此时回车键作为分行标志， A=[1,2,3 4,5,6 7,8,9] （3）使用 M 文件创建大矩阵，当矩阵维数非常大时，可以创建 m 文件，在 m 文 件中输入数据或者导出数据文件。 2、矩阵的基本运算 ①矩阵的四则运算。其中乘法运算要注意相乘的双方有相邻公共维，除法分为左除 “\” (A\B=inv(A)*B)和右除”/”（ A/B=A*inv (B)）(需要计算逆矩阵) ②矩阵的逆运算。inv 函数。 ③矩阵的幂运算。^。 ④矩阵的指数运算。exp(返回每个元素的指数值)，expm([V,D] = EIG(X) and EXPM(X) = V*diag(exp(diag(D)))/V)，expm1(exp(x)-1) ⑤矩阵的对数运算。logm ⑥矩阵的特征值函数。eig 和 eigs(适合于大型稀疏方阵) ⑦矩阵的奇异值函数。svd（[U,S,V] = SVD(X)，X = U*S*V'）和 svds ⑧矩阵的条件数函数。cond（矩阵 A 的条件数等于 A 的范数与 A 的逆的范数的 乘积，c = cond(A,p)等价于 norm(A,p) * norm(inv(A),p)） ，condest(1 范数的条件数的估计 值)， rcond ⑨特征值的条件数函数。codeig（[V,D,s] = condeig(A) 等价于[V,D] = eig(A); s6 =condeig(A);） ⑩范数函数。norm（1-范数：即列范数，矩阵的各列绝对值之和的最大值；2-范数： 所有元素的平方和开根号(默认)；无穷范数：即行范数，矩阵各行的绝对值之和的最大 值） ， normest (矩阵的 2 范数的估计值) 其他还有秩函数 rank， 迹函数 trace， 零空间函数 null （又称为核空间， X=null(A)， 则 A*X=0，X?*X=I） ，正交空间函数 orth(B = orth(A),则 B'*B = eye(rank(A)))，伪逆函数 pinv 等。 3、特殊向量和特殊矩阵 （1）特殊向量 t=[0:0.1:10] %产生从 0 到 10 的行向量，元素之间间隔为 0.1 t=linspace(n1,n2,n) %产生 n1 和 n2 之间线性均匀分布的 n 个数 (缺省 n 时,产生 100 个数) t=logspace(n1,n2,n) (缺省 n 时,产生 50 个数) %在和之间按照对数距离等间距产生 n 个数。 （2）特殊矩阵 i)单位矩阵 eye(m), eye(m,n) 可得到一个可允许的最大单位矩阵而其余处补 0, eye(size(a)) 可以得到与矩阵 a 同样大小的单位矩阵。 ii) 所有元素为 1 的矩阵 ones(n)，ones(size(a))，ones(m， n)。 iii）所有元素为 0 的矩阵 zeros(n)，zeros(m,n)。 iv）空矩阵是 一个特殊矩阵，这在线性代数中是不存在的。例如 q=[ ] 矩阵 q 在工作空间之中，但它的大小为零。通过空矩阵的办法可以删除矩阵的行与 列。例如 a(:,3)=[] 表示删除矩阵 a 的第 3 列。 v）随机数矩阵7 rand(m,n) 产生 m× n 矩阵，其中的元素是服从[0,1]上均匀分布的随机数。 randint(m,n,[min,max]) 产生 m× n 矩阵，其中的元素是[min,max]上的随机整数。 normrnd(mu,sigma,m,n)产生 m× n 矩阵，其中的元素是服从均值为 mu，标准差 为 sigma 的正态分布的随机数。 exprnd(mu,m,n) 产生 m× n 矩阵，其中的元素是服从均值为 mu 的指数分布的随 机 数。 poissrnd(mu,m,n) 产生 m× n 矩阵，其中的元素是服从均值为 mu 的泊松 （Poisson） 分布的随机数。 unifrnd(a,b,m,n) 产生 m× n 矩阵， 其中的元素是服从区间[a,b]上均匀分布的随机数。 r = mvnrnd(MU,SIGMA,cases) 产生 cases 对均值向量为 MU，协方差阵为 SIGMA 的多维正态分布的随机数。 vi）随机置换 randperm(n)产生 1 到 n 的一个随机全排列。 perms([1:n])产生 1 到 n 的所有全排列。 vii) 稀疏矩阵 稀疏矩阵是指矩阵中零元素很多，非零元素很少的矩阵。对于稀疏矩阵，只要存放 非零元素的行标、列标、非零元素的值即可，可以按如下方式存储 （非零元素的行地址，非零元素的列地址），非零元素的值。 在 Matlab 中无向图 和有向图邻接矩阵的使用上有很大差异。 对于有向图，只要写出邻接矩阵，直接 使用 Matlab 的命令 sparse 命令，就可以把 邻接矩阵转化为稀疏矩阵的表示方式。 对于无向图，由于邻接矩阵是对称阵，Matlab 中只需使用邻接矩阵的下三角元素， 即 Matlab 只存储邻接矩阵下三角元素中的非零元素。 稀疏矩阵只是一种存储格式。Matlab 中，普通矩阵使用 sparse 命令变成稀疏矩阵， 稀疏矩阵使用 full 命令变成普通矩阵。 例 1 a=zeros(5); a(1,[2,4])=[3,4]; a(3,[2:4])=[1 3 8]; a(5,[1,5])=[6,7]8 b=sparse(a) %普通矩阵转化成稀疏矩阵 c=full(b) %稀疏矩阵转化成普通矩阵 其他一些特殊矩阵如下表所示： 函数 compan gallery hadamard hankel hilb invhilb 功能 伴随阵 Higham 测试阵 Hardamard 矩阵 Hankel 矩阵 Hilbert 矩阵 反 Hilbert 矩阵 函数 magic rosser toeplitz pascal vander wilkinson 功能 魔方阵 经典对称特征值测试阵 Toeplitz 矩阵 Pascal 矩阵 范德蒙矩阵 Wilkinson’s 特征值测试矩阵4、矩阵的特殊操作 ①变维。 有两种方法，使用冒号（:）和使用函数 reshape 使用“:”表达式对两个矩阵进行变维操作，需要预先定义两个矩阵的维数（例如 A(:,3)=[]） ；reshape 有两种形式，分别为 reshape(X,M,N)和 reshape(X,M,N,P…) ②变向 主要函数如下表所示： 函数 fiplr fipud fipdim Rot90 ③矩阵的抽取 对角线元素抽取函数 diag(X,k)/diag(v,k)，抽取矩阵 X 的第 k 条对角线的元素向量/ 使得向量 v 为所得矩阵的第 k 条对角线元素。 上三角元 素抽取 tril(X,k)和下三角元素抽取 triu(X,k) ④扩展 两种方法：利用对矩阵标示块的赋值命令 X(m1:m2,n1:n2)=a 生成大矩阵， 其中 m2-m1+1 必须等于 a 的行维数，n2-n1+1 必须等于 a 的列维数，生成 m2 × n2 维的 功能 矩阵左右翻转 矩阵上下翻转 矩阵特定维翻转 矩阵反时针 90 翻转 函数 diag tril triu 功能 产生或提取对角阵 产生下三角 产生上三角9 矩阵 X；利用小矩阵组合生成大矩阵，要严格注意矩阵大小的匹配。 （示例： X(1:10,1:10)=zeros(10,10)，LX=[X,X;X,X]） Matlab 中冒号（： ）的使用方法小结： （1）用于生成向量，a:b，一般要求 a&b，否则生成空矩阵。 （2）如果 b-a 不是整数时，则向量的最后一位数为 n+a，其中 n=fix(b-a)(向零取整) （3）a:c:b 表示[a,a+c,…,a+n*c]，其中 n=fix((b-a)/c)，当 c&0 且 a&b 或者 c&0 且 a&b 时出现空矩阵。 （4）A(:)以一列的方式显示 A 中所有的元素 （5）b=A(i,:)表示将 A 中第 i 行存入 b 中 （6）b=A(:,j)表示将 A 中第 j 列存入 b 中 （7）b=A(j:k)表示将 A 中第 j 到第 k 个元素存入 b 中（Matlab 中矩阵按列存储） （8）b=A(:,c:d)表示将 A 中第 c 列到第 d 列存入 b 中，要求 c，d 不能超过 A 的列 数。 （9）当矩阵很大时，不知道矩阵的维数，可以使用 end 作为矩阵的最后一行或者 一列或者最后一个元素。例如 b=A(1:2,2:end)获取矩阵 A 右上角的元素。2.4 数组及其运算数组与矩阵在形式上完全一致，只是运算与矩阵不同。同型矩阵之间的运算通常称 为数组运算。 （矩阵的数组运算） 1、基本数组运算 ①四则运算。数组的乘除法是指两个同维数组间对应元素之间的乘除法，运算符 为”.*”,”./”和”.\”。数组与常数之间的运算可以加”.”，也可以不加。 ②幂运算。.^,对每个数组元素的幂运算。 ③指数运算 exp,对数运算 log 和开方运算 sqrt。 2、数组函数运算 只要把运算的数组带入到函数中就可以了，通用 形式为 funname(A) 3、数组的逻辑运算和关系运算 指令 & 含义 小于 函数名 lt10 &= & &= == ~= & | ~小于等于 大于 大于等于 等于 不等于 逻辑 与 逻辑 或 逻辑 非le gt ge eq ne and or not指令 xor any all isnan含义 不相同就取 1，否则取 0 只要有非 0 就取 1，否则取 0 全为 1 取 1，否则为 0 为数 NaN 取 1，否则为 0指令 isequal ismember isempty isletter含义 相等取 1，否则取 0 两个矩阵是属于关系取 1， 否则取 0 矩阵为空取 1，否则取 0 是字母取 1， 否则取 0 （可以是字符 串）isinf isfinite ischar find isnumeric为数 inf 取 1，否则为 0 有限大小元素取 1，否则为 0 是字符串取 1，否则为 0 寻找非零元素坐标 判断数值矩阵isstudent isprime isreal isspace islogical学生版取 1 质数取 1，否则取 0 实数取 1，否则取 0 空格位置取 1，否则取 0 判断逻辑数组11 3、Matlab 中的 M 文件及程序调试M 文件分为两种，一种是脚本文件，由一系列 Matlab 的命令组成，可以直接运行； 一种是函数文件，必须由其他 M 文件或者在命令行窗口中调用执行。函数文件具有一 定的通用性，并且可以进行递归调用。3.1、Matlab 中流程控制语句（1）if 语句 有三种形式：if(表达式) 语句组 A; end if(表达式) 语句组 A，else 语句组 B，end if(表达式 1) 语句组 A，elseif(表达式 2) 语句组 B, else 语句组 C，end （2）循环语句 while（表达式）语句组 A , end for k=初值：增量：终值 语句组 A，end 循环语句可以结合 break 命令来控制程序的执行顺序。 （3）switch 语句 switch 表达式(标量或者字符串) case 值 1 语句组 A case 值 2 语句组 B … otherwise 语句组 N end 用法示例： 示例一：输入数 n，判断其奇偶性 解 Matlab 程序如下： n=input(?n=?); if rem(n,2)==0 %求余数，与 mod 函数的区别，当同符号时两者相同，否则不同 A=?even?; else A=?odd?,12 end 示例二：求 Matlab 中的最大实数 解：Matlab 程序如下： x=1; while x~=inf, x1=x; x=2*x; %为了获取更接近 realmax 的值，可以改为 x=1.01*x end x1 示例三：求 Matlab 的相对精度 解：matlab 程序如下： y=1; while 1+y&1 y1=y; y=y/2; end y13.2 脚本文件主要特征如下： (1)一般由 clc,clear 命令开始，清除掉屏幕和工作空间中原有的变量和图形，以避免 其他已执行的程序残留数据对本程序的影响。程序中应该添加有注释。 (2)接下来是程序的主体，如果文件中有全局变量，则需要使用 Global 在程序的起 始部分注明。为了提高程序的可读性，注意语句的缩进。 (3)整个程序应按照 Matlab 标识符的要求起文件名，扩展名为 m。 示例：列出求素数的程序 解 Matlab 程序如下： clear,clc N=input('N='),x=2:N; %列出从 2 到 N 的全部自然数 for u=2:sqrt(N) %依次列出除数 n=find(rem(x,u)==0&x~=u);%找到能被 u 整除而不等于 u 的数序号 x(n)=[]; %剔除该数 end x13 3.3 函数文件函数文件与脚本文件的区别： (1)由 function 开头，后跟的函数名必须与文件名相同； (2)有输入变量和输出变量，可进行变量传递； (3)除非使用 Global 声明，程序中的变量均为局部变量，运行后不保存在工作空间 中。 示例 1：使用函数文件写出求素数的程序 解：Matlab 程序如下： function y=qiuprime(N) %求出 1 到 N 之间的素数，返回素数组成的向量 x=2:N; %列出从 2 到 N 的全部自然数 for u=2:sqrt(N) %依次列出除数 n=find(rem(x,u)==0&x~=u);%找到能被 u 整除而不等于 u 的数序号 x(n)=[]; %剔除该数 end y=x; 示例 2：写一个递归函数求 n! 解： function f=factor(n) if n&=1 f=1; else f=factor(n-1)*n; end end3.4 函数参数的可调性Matlab 中有两个永久变量 nargin 和 nargout 分别记录调用该函数时的输入实参和输 出实参的个数。只要在函数中包含这两个变量，就可以准确的知道该函数文件被调用时 的输入输出参数个数，从而决定函数如何进行处理。其他类似的变量还有 varargout， varargin。14 示例 1： function fout=examp(a,b,c) if nargin==1 fout=a; elseif nargin==2 fout=a+b; elseif nargin==3 fout=a+b+c; end 示例 2： function [x,y,z]=examp(a,b,c) if nargout==1 x=a; elseif nargout==2 x=a; y=b; elseif nargout==3 x=a; y=b; z=c; end end 示例 3： function vartest(argA, argB, varargin) optargin = size(varargin,2); %可选输入 stdargin = nargin - %标准输入 fprintf('Number of inputs = %d\n', nargin) fprintf(' Inputs from individual arguments(%d):\n', ... stdargin) if stdargin &= 1 fprintf(' %d\n', argA) end if stdargin == 2 fprintf(' %d\n', argB) end fprintf(' Inputs packaged in varargin(%d):\n', optargin) for k= 1 : size(varargin,2) fprintf(' %d\n', varargin{k}) end end15 3.5 程序调试M 文件编辑器中提供了强大的程序调试功能，使用方式如同 VC 程序中的调试器。 在 M 文件编辑器中的工具栏中有如下程序调试图标：这些图标的使用和菜单 Debug 中的一些命令相同。16 4、Matlab 绘图命令4.1 二维绘图命二维绘图的基本命令有 plot，loglog，semilogx，semilogy 和 polar。它们的使用方法 基本相同，其不同特点是在不同的坐标中绘制图形。plot 命令使用线性坐标空间绘制图 形；loglog 命令在两个对数坐标空间中绘制图形；而 semilogx(或 semilogy)命令使用 x 轴（或 y 轴）为对数刻度，另外一个轴为线性刻度的坐标空间绘制图形；polar 使用极 坐标空间绘制图形。 二维绘图命令 plot 为了适应各种绘图需要，提供了用于控制线色、数据点和线型的 3 组基本参数。它的使用格式如下： plot(x,y,‘color_point_linestyle’) 该命令是绘制 y 对应 x 的轨迹的命令。y 与 x 均为向量，且具有相同的元素个数。 用字符串，color_point_linestyle? 完成对上面 3 个参数的设置。 线色：r-red，g-green，b-blue，w-white，k-black，i-invisible，y-yellow，c-cyan （青色） ，m-紫色。 数据点：.（圆点） ，o（小圆圈） ，x（叉号） ，+（加号） ，*（星号） ，s （square 方形） ， h（hexagram 六角星） ，d（diamond 菱型） ，p（pentagram 五角星） ，v（下三角） ，^（上 三角） ，&（左三角） ，&（右三角） 。 线型：-（实线） ，--（虚线） ，-.（点画线） ，:（点线） 。 当 plot(x,y)中的 x 和 y 均为 m× n 矩阵时，plot 命令将绘制 n 条曲线。 plot(t,[x1,x2,x3])在同一坐标轴内同时绘制三条曲线。 如果多重曲线 对应不同的 x 轴向量绘制，可使用命令 plot(t1,x1,t2,x2,t3,x3) 式中 x1 对应 t1，x2 对应 t2 等等。在这种情况下，t1，t2 和 t3 可以具有不同的元素 个数，但要求 x1，x2 和 x3 必须分别与 t1，t2 和 t3 具有相同的元素个数。 subplot 命令使得在一个屏幕上可以分开显示 n 个不同坐标系，且可分别在每一个 坐标系中绘制曲线。其命令格式如下 subplot(r,c,p)17 该命令将屏幕分成 r× c 个子窗口，而 p 表示激活第 p 个子窗口。窗口的排号是从 左 到右，自上而下。 下面对几个特殊的坐标系统 进行简要介绍： ①对数坐标曲线，主要有 semilogx,semilogy 和 loglog，前两个分别以 x 坐标和 y 坐标为对数坐标，后一个是双对数坐标。 例如： x=1:0.1*pi:2* y=sin(x); semilogx(x,y,'-*') x = 0:.1:10; semilogy(x,10.^x) x = logspace(-1,2); loglog(x,exp(x),'-s') grid on ②双纵坐标（双 y 轴坐标系）函数 plotyy，调用形式为： plotyy(X1,Y1,X2,Y2) plotyy(X1,Y1,X2,Y2,fun) fun 可以是 plot、semilogx、semilogy 或 loglog plotyy(X1,Y1,X2,Y2,fun1,fun2) fun1 绘制(X1,Y1)，fun2 绘制(X2,Y2) ③极坐标绘图函数 polar 调用形式为：polar(theta,rho) 例 1：绘制极坐标下的平面曲线 ? ? a cos(b ? n? ) ，并讨论参数 a,b,n 对曲线的影响。 解 Matlab 程序如下： N=100; theta=linspace(0,2*pi,N); for i=1:2 a(i)=input('a='); b(i)=input('b='); n(i)=input('n='); rho(i,:)=a(i)*cos(b(i)+n(i)*theta); subplot(1,2,i),polar(theta,rho(i,:)); end 几个比较漂亮的极坐标系下的图 形： (1)蝴蝶图案18 t=0:0.01:36; f=exp(cos(t-pi/2))-2*cos(4*(t-pi/2))+sin((t-pi/2)/12).^5; polar(t,f,'r') (2)枫叶图案 t=-pi/2:0.05:1.5* f=100./(100+(t-pi/2).^8)*2-sin(7*t)-cos(30*t)/2; polar(t,f,'r') 在图形绘制完毕后， 执行如 下命令可以再在图中加入标题、 标号、 说明和分格线等。 这些命令有 title，xlabel，ylabel，text，gtext，legend 等。它们的命令格式如下 title(“My Title”)，xlabel(“My X-axis Label ”)，ylabel(“My Y-axis Label”)， xlabel({'first line';'second line'})， text(x, y,'Text for annotation')， gtext('Text for annotation')， grid on gtext 命令是使用鼠标器定位的文字注释命令。当输入命令后，可以在屏幕上得到 一个光标， 然后使用鼠标器控制它的位置。 按鼠标器的左键， 即可确定文字设定的位置。 hold on 是图形保持命令，可以把当前图形保持在屏幕上不变，同时在这个坐标系 内绘制另外一个图形。 hold on 命令是一个交替转换命令， 即执行一次， 转变一个状态 （相 当于 hold on、hold off） 。 还可以设置坐标轴的范围，使用命令 axis，其格式为： axis([xmin xmax ymin ymax]) axis([xmin xmax ymin ymax zmin zmax cmin cmax]) 示例： x = 0:.01:pi/2; plot(x,tan(x),'-ro') 图形效果如下所示：默认情况下，横纵坐标轴的范围是根据函数自变量和因变量的值自动变化的，有时 效果不好，此时需要设定横纵坐标轴的范围：axis([0 pi/2 0 5])，效果如上右图所示。19 解 画图的 Matlab 程序如下，画出的图见图 1。 clc, clear x=-2*pi:0.1:2* y1=sin(x); y2=sin(x+pi/3)+2; y3=cos(x); plot(x,y1,'.-'); hold on %图形保持命令 plot(x,y2,'*-'); plot(x,y3,'-o'); h=legend('sin($x$)','sin($x+\frac{\pi}{3}$)','cos($x$)') %latex 格式显示 set(h,'Interpreter','latex') %设置 Interpreter 的属性值为 latex，可以使用数学公式 xlabel('$x$','Interpreter','latex') %latex 格式显示 ylabel('$y$','Interpreter','latex') %latex 格式显示图1解 画图的 Matlab 程序如下，画出的图形结果见图 2。 clc, clear x=-2*pi:0.1:2* y1=sin(x); y2=sin(x+pi/3)+2; y3=cos(x);20 subplot(3,1,1), plot(x,y1,'.-'), title('sin($x$)','Interpreter','latex') subplot(312), plot(x,y2,'*-') title('sin($x+\frac{\pi}{3}$)','Interpreter','latex') ylabel('$y$','Interpreter','latex') %latex 格式显示 subplot(313), plot(x,y3,'-o') title('cos($x$)','Interpreter','latex') %latex 格式显示 xlabel('$x$','Interpreter','latex') %latex 格式显示4.2 复数绘图plot(z) ：z 为复数时相当于 plot(real(z),imag(z))；如果是双变量如 plot(t,z)，则 z 中 的虚数部分将丢弃。在复平面中绘出多条曲线，必须使用 hold on 命令，或者把多条曲 线的实部和虚部明确的写出。即 plot(real(z1),imag(z1),real(z2),imag(z2)) Matlab 中专门用来绘制复变量函数图形的相关命令是 cplxmap， cplxgrid， cplxroot， 格式如下： z = cplxgrid(m) %产生(m+1)*(2m+1)的极坐标下的复数数据网格 cplxmap(z,f(z),(optional bound)) %绘制复变函数的图形，以 xy 平面为自变量所在的 复平面，以 z 轴表示复变函数的实部，颜色表示复变函数的虚部。 cplxroot：画复数的 n 次函数曲面21 cplxroot(n) %画复数 n 次根的函数曲面,复数为最大半径为 1 的圆面 cplxroot(n,m) %画复数 n 次根的函数曲面,复数为最大半径为 1 的圆面,为（m+1）* （2m+1）的方阵 例 1：绘制 z=exp((-0.1+i)*t)的复数图形 解 Matlab 程序如下： t=0:0.1:15; z=exp((-0.1+i)*t); subplot(2,2,1) plot(z), % pause 可以暂停程序的执行， 通过按任意键可使程序继续进行， 也 可以在 pause(n)中设置时间使执行结果出现动态效果。 title('复数绘图 plot(z)'); subplot(2,2,2) plot(t,z),pause title('plot(t,z)') subplot(2,2,3) polar(angle(z),abs(z)); title('polar(angle(z),abs(z))') subplot(2,2,4) semilogx(t,z); title('semilogx(t,z)') 注：Matlab 中使用函数 C=complex(A,B)构造复数。 例 2：绘图 y=cos(x+i)的图形 解法一：x=-pi:0.1: fun=@(x)cos(x+i); %如果直接使用 plot(x,y)，会忽略虚部。 plot(fun(x)) 解法二：plot(cos(x+i)); 例 3：绘制 f=z^4 的图形，其中 z 为复数 解：z=cplxgrid(30); cplxmap(z,z.^4)4.3 显函数，符号函数或隐函数的绘图fplot(fun,lims)绘制由字符串 fun 指定函数名的函数在 x 轴区间为 lims=[xmin, xmax] 的函数图。若 lims=[xmin,xmax,ymin,ymax]，则 y 轴也被限制。22 解 （1）首先用 M 文件 fun1.m 定义函数 f(x)如下 function y=Afun1(x); if x&1 y=x+1; else y=1+1./x; end 在 matlab 命令窗口输入 fplot('Afun1',[-3,3]) 就可画出 函数 f (x) 的图形。 （2）可以使用匿名函数，编写程序如下 fun2=@(x) (x+1)*(x&1)+(1+1/x)*(x&=1); fplot(fun2,[-3,3]) ezplot(f)绘制符号函数或者隐函数 f(x)的图形，x 轴的近似范围为[2?,?2??。 ezplot(f,[xmin,xmax])使用输入参数来代替默认横坐标范围[2?,?2?? 。 ezplot 函数的其 他格式有： ? ezplot(fun2)绘制 fun2(x,y)=0 的隐函数曲线，默认 x,y 的范围是[2 ?,?2??? ezplot(fun2,[xmin,xmax,ymin,ymax]) 示例： function z = myfun(x,y,k) %建立 M 文件 z = x.^k - y.^k - 1; %在命令窗口中输入：ezplot(@(x,y)myfun(x,y,2));解 ezplot('cot(x)')解 ezplot('x^2+y^2/4=1') 例 6-1：同一坐标系下绘制曲线 x^2+y^2=1 和 x^2-y^2=1 的曲线 解：ezplot('x^2-y^2=1');23 ezplot('x^2+y^2=1'); colormap([0,0,1]);4.4 三维图形在实际工程计算中，最常用的三维绘图是三维曲线图、三维网格图和三维曲面图 3 种基本类型。 与此对应， Matlab 也提供了一些三维基本绘图命令， 如三维曲线命令 plot3， 三维网格图命令 mesh 和三维表面图命令 surf。 1. 三维曲线 plot3(x,y,z)通过描点连线画出曲线， 这里 x,y,z 都是 n 维向量， 分别表示该曲线上点 集的横坐标、纵坐标、竖坐标。 例7 t=0:0.05:20* x=sin(t); y=cos(t); z=t.*sin(t).*cos(t); plot3(x,y,z), title('Line in 3-D Space') xlabel('X'), ylabel('Y'), zlabel('Z'), grid on 2. 三维网格图 命令 mesh(x,y,z)画网格曲面。这里 x,y,z 是三个同维数的数据矩阵，分别表示数据 点的横坐标、纵坐标、竖坐标，命令 mesh(x,y,z)将该数据点在空间中描出，并连成网格。解 x=-3:0.1:3;y=-5:0.1:5; x1=ones(size(y'))*x;y1=y'*ones(size(x)); [x2,y2]=meshgrid(x,y); %meshgrid 函数生成 2D/3D 网格矩阵，用于分割空间，在绘 制 3D 网格图或者表面图时都需要调用该函数对生成绘图时所用的数据。 z1=(sin(x1.*y1)+eps)./(x1.*y1+eps); z2=(sin(x2.*y2)+eps)./(x2.*y2+eps);24 subplot(1,2,1),mesh(x1,y1,z1) subplot(1,2,2),mesh(x2,y2,z2) 例 8-1：绘制两个空间相交的平面 z=x+2y-1 和 z=x-y+2 解 ： x=-3:0.1:3;y=-3:0.1:3; [X,Y]=meshgrid(x,y); mesh(X,Y,X+2*Y-1) hold on mesh(X,Y,X-Y+2) 3. 表面图 命令 surf(x,y,z)画三维表面图，这里 x,y,z 是三个同维数的数据矩阵，分别表示数据 点的横坐标、纵坐标、竖坐标。解 [x,y]=meshgrid([-3:0.2:3]); z=(sin(x.*y)+eps)./(x.*y+eps); surf(x,y,z) 4. 旋转曲面 方法一：使用命令函数 cylinder [X,Y,Z] = cylinder(r) 这里的 r 表示构成旋转曲面的曲线。解 Matlab 程序如下。 x=0:10:600; [X,Y,Z]=cylinder(30*exp(-x/400).*sin((x+25*pi)/100)+130); surf(X,Y,Z) 方法二：将旋转曲面用参数方程表示。25 解 因为这里的函数是隐函数，化成显函数后有两支，必须使用参数方程，旋转面 的参数方程为画图的 Matlab 程序如下： alpha=[0:0.1:2*pi]'; beta=0:0.1:2* x=4*cos(alpha)*ones(size(beta)); y=(5+4*sin(alpha))*cos(beta); z=(5+4*sin(alpha))*sin(beta); surf(x,y,z) 或者利用绘制三维隐函数曲面图形的命令 ezsurf 或者 ezmesh，其命令格式如下： ezmesh(fun) ezmesh(fun,domain) domain 要求是一个向量，指定 x,y 轴的范围 ezmesh(funx,funy,funz) %参数方程绘图 funx(s,t), funy(s,t), and funz(s,t) over the square: -2π & s & 2π, -2π & t & 2π. ezsurf(fun) ezsurf(fun,domain) ezsurf(funx,funy,funz) 画图的 Matlab 程序也可以写成 x=@(alpha,beta) 4*cos(alpha); y=@(alpha,beta) (5+4*sin(alpha))*cos(beta); z=@(alpha,beta) (5+4*sin(alpha))*sin(beta); ezsurf(x,y,z) 对于形如 f(x,y,z)=0 的三维隐函数的图像没有现成的函数可以画，但可以利用 isosurface 函数绘制三角网格图。例如 f=@(x,y,z)x.^2+y.^2+z.^2-10;%定义函数 f=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,0);%用 isosurface 得到函数 f=0 图形的点和面26 patch('faces',p,'vertices',v,'facevertexcdata',jet(size(v,1)),'facecolor','w','edgecolor','flat'); %用 patch 绘制三角网格图并设定色彩 view(3); axis equal 5、绘制柱面 柱面平行于某个坐标轴，方程中不出现某个坐标轴的变量，方程表示为 F(x,y)=0 或 者 F(x,z)=0 或者 F(y,z)=0 示例一：画出方程 z ? 2 ? x2 表示的柱面 解：方程中的 x 是自变量矩阵，z 是因变量，则另一个自变量矩阵为 y，自变量平 面是 xoy 面，x 轴是真正的自变量，y 轴是柱面方向。Matlab 程序如下： u=linspace(-5,5,10)';%设定参数列向量 u v=linspace(-5,10,10); %设定参数行向量 v X=u*ones(size(v)); %构成自变量矩阵 X Y=ones(size(u))*v; %构成自变量矩阵 Y Z=2-X.^2; %求因变量 Z mesh(X,Y,Z) 示例二：画出方程 ? x 2 ? y2 ? 1 表示的柱面 4解：方程整理为显示函数： y ? ?2 1 ? x2 ，分别绘制正负两个分量。 clear,clc u=linspace(-5,5,10)';%设定参数列向量 u v=linspace(-5,10,10); %设定参数行向量 v X1=u*ones(size(v)); %构成自变量矩阵 X1 Y1=2*sqrt(1+X1.^2); %求因变量 Y1 Y2=-2*sqrt(1+X1.^2); Z1=ones(size(u))*v; %求因变量 Z1 mesh(X1,Y1,Z1),hold on mesh(X1,Y2,Z1) 6. 其它二次曲面 对于旋转面，如果母线的方程可以表示成关于旋转轴变量的显式函数，则可以直接 使用 Matlab 工具箱中的命令 cylinder，否则必须把旋转面化成参数方程，然后使用 ezmesh 或 ezsurf 命令绘图。对于其它的二次曲面，如果可以写成显函数直接使用命令 ezmesh 或 ezsurf，否则必须先化成参数方程。还有一些特殊的绘制函数如27 cylinder,ellipsoid 等命令。解： （1）x=@(s,t) 3*sec(s)*cos(t); %化为参数方程 y=@(s,t) 3*sec(s)*sin(t); z=@(s,t) 2*tan(s); ezmesh(x,y,z) （2）x=@(s,t) 3*sec(s); y=@(s,t) 2*tan(s)*cos(t); z=@(s,t) 2*tan(s)*sin(t); ezmesh(x,y,z) （3）ezsurf(@(y,z) y.^2,50) %直接调用 ezsurf （4）x=@(s,t) 3*tan(s)*cos(t); y=@(s,t) 2*tan(s)*sin(t); z=@(s,t) tan(s); ezsurf(x,y,z) （5）ellipsoid(0,0,0,3,2,sqrt(6)) %专门绘制椭球面 （6）ezsurf(@(x,y)x*y) （7）x=@(s,t) 3*cos(s);28 y=@(s,t) 2*sin(s); z=@(s,t) ezmesh(x,y,z) 7、绘制空间两曲面的交线 示例：绘制由水平截面与方程 z ? x2 ? 2 y 2 构成的马鞍面形成的交线，并讨论等高线 和方向导数（梯度）的意义 解： 水平平面与曲面的交线就是等高线， 在 Matlab 中绘制等高线有两个命令， contour 和 contour3，前者把等高线画在 xoy 平面上，后者把等高线画在一定高度的平面上，使 之成为立体的，与所在曲面对应。 Matlab 程序如下： clc,clear [x,y]=meshgrid(-10:.2:10);%确定计算和绘图的定义域网格 z1=(x.^2-2*y.^2)+ %第一个曲面方程 a=input('a=(-50&a&50)'); z2=a*ones(size(x)); %水平面方程 z2=a，z2 必须与 x，y 具有相同的维数， subplot(1,3,1),mesh(x,y,z1);mesh(x,y,z2);%分别画出两个曲面 v=[-10 10 -10 10 -100 100];axis(v), %确定第一个分图的坐标系 colormap(gray); r0=abs(z1-z2)&=1; %求两曲面 z 坐标只差小于 0.5 的网格 zz=r0.*z2;yy=r0.*y;xx=r0.*x; %求这些网格上的坐标值，即交线坐标值 subplot(1,3,2),plot3(xx,yy,zz,'x'); %画出这些点 axis(v), %使第二个分图取第一个分图的坐标系 pause,subplot(1,3,3); contour3(x,y,z1,20); %用等高线命令求出 20 条不同高度的交线。 等高线与方向导数和梯度的概念密切相关，函数 z1 在每一点的梯度与该处的等高 线垂直，也就是指向最陡的方向。Matlab 中求梯度用 gradient 函数，quiver 函数画出梯 度向量，这两个函数要求在给定的点阵上求梯度和画梯度向量。Matlab 程序如下： [x,y]=meshgrid(-10:2:10); z1=(x.^2-2*y.^2)+ contour(x,y,z1,20); [px,py]=gradient(z1,2,2); %以步长为 2 求 z1 的梯度的 x，y 分量 quiver(x,y,px,py); %绘制梯度向量 8、离散点的图形绘制相关函数 griddata，TripScatterInterp，scatter（scatter3） griddata 用来对离散数据进行曲面拟合 (1)ZI = griddata(x,y,z,XI,YI) 用二元函数 z=f(x,y)的曲面拟合有不规则的数据向量29 x,y,z。griddata 将返回曲面 z 在点（XI,YI）处的插值。曲面总是经过这些数据点（x,y,z） 的。输入参量（XI,YI）通常是规则的格点（像用命令 meshgrid 生成的一样） 。XI 可以 是一行向量，这时 XI 指定一有常数列向量的矩阵。类似地，YI 可以是一列向量，它 指定一有常数行向量的矩阵。 (2)[XI,YI,ZI] = griddata(x,y,z,xi,yi) 返回的矩阵 ZI 含义同上， 同时， 返回的矩阵 XI,YI 是由行向量 xi 与列向量 yi 用命令 meshgrid 生成的。 (3)[XI,YI,ZI] = griddata(.......,method) 用指定的算法 method 计算： ?linear ?：基于三角形的线性插值（缺省算法） ； ?cubic?： 基于三角形的三次插值； ?nearest?：最邻近插值法； 示例： x = rand(100,1)*4-2; y = rand(100,1)*4-2; z = x.*exp(-x.^2-y.^2); ti = -2:.25:2; [xi,yi] = meshgrid(ti,ti); zi = griddata(x,y,z,xi,yi); mesh(xi,yi,zi), hold on, plot3(x,y,z,'o'), hold off 在新版本的 Matlab 中，该函数逐渐被 TriScatteredInterp 函数所代替， TriScatteredInterp 的示例如下： x = rand(100,1)*4-2; y = rand(100,1)*4-2; z = x.*exp(-x.^2-y.^2); F = TriScatteredInterp(x,y,z); ti = -2:.25:2; [qx,qy] = meshgrid(ti,ti); qz = F(qx,qy);%计算指定位置的插值 mesh(qx,qy,qz); plot3(x,y,z,'o'); Matlab 中绘制散点图的命令是 scatter（scatter3） ，命令格式如下： scatter(X,Y,S,C) %X，Y 对应散点值，S 代表标记大小，C 代表颜色值，它都可以是 向量。 scatter3(X,Y,Z,S,C) 示例：30 A=[1.486,3.059,0.1;2.121,4.041,0.1;2.570,3.959,0.1;3.439,4.396,0.1; 4.505,3.012,0.1;3.402,1.604,0.1;2.570,2.065,0.1;2.150,1.970,0.1; 1.794,3.059,0.2;2.121,3.615,0.2;2.570,3.473,0.2;3.421,4.160,0.2; 4.271,3.036,0.2;3.411,1.876,0.2;2.561,2.562,0.2;2.179,2.420,0.2; 2.757,3.024,0.3;3.439,3.970,0.3;4.084,3.036,0.3;3.402,2.077,0.3; 2.879,3.036,0.4;3.421,3.793,0.4;3.953,3.036,0.4;3.402,2.219,0.4; 3.000,3.047,0.5;3.430,3.639,0.5;3.822,3.012,0.5;3.411,2.385,0.5; 3.103,3.012,0.6;3.430,3.462,0.6;3.710,3.036,0.6;3.402,2.562,0.6; 3.224,3.047,0.7;3.411,3.260,0.7;3.542,3.024,0.7;3.393,2.763,0.7]; x=A(:,1);y=A(:,2);z=A(:,3); scatter(x,y,5,z)%散点图4.5 动态可视化图形Matlab 中的动画命令，moviein，getframe 和 movie，用 getframe 把 Matlab 产生的 图形存储下来，每个图形成一个很大的列向量；再用 N 行这样的列保存 N 幅图，成为 一个大矩阵；movie 命令把他们连接起来重放，产生动画效果；moviein 用来预留存储空 间以加快运行的速度。有些情况也需要借助 pause 命令和循环语句来实现点轨迹的动态 演示。绘制彗星图的命令 comet 或者 comet3 也可以实现点的运动轨迹，不过这种方式 速度较快。 示例 1： axis equal, %把坐标设成相等比例 M=moviein(16);%为变量 M 预留 16 幅图的存储空间 for j=1:16 plot(fft(eye(j+16))); M(:,j)= end movie(M,30) %以每秒 30 帧的速度播放 M 中的图形 示例 2：使用 comet/comet3 绘制彗星图 （1） t = 0:.01:2* %返回当前坐标轴下的一帧图像31 x = cos(2*t).*(cos(t).^2); y = sin(2*t).*(sin(t).^2); comet(x,y); （2） n=10; t=n*pi*(0:0.0005:1); x=sin(t); y=cos(t); plot(x,y,'g'); comet(x,y,0.01); 示例 3：极坐标下一个点沿着图形移动的轨迹 a=0:.01:2* b=3; polar(a,b*(1-cos(a)),'*r'); for a=0:0.05:2*pi cla(findobj(gca,'color','r')); %findobj 函数用来使用特定的属性获取图形句柄， cla 用来清理由函数句柄指定的图形所在的坐标轴。 h2=polar(a,b*(1-cos(a)),'b'); set(h2,'linestyle','.','markersize',30); %设置句柄指定函数图形的属性 pause(0.001) end 示例 4： 普通二维坐标下一个点沿着曲线移动的轨迹， 程序与上面的结果基本相同。 x=-2*pi:0.1:2* y=sin(x); plot(x,y,'b-'); hold on for m=-2*pi:0.05:2*pi cla(findobj(gca,'color','b')); %cla 清除当前坐标轴上面的子对象 h=plot(m,sin(m),'r'); set(h,'LineStyle','.','Marker','*','MarkerSize',20); pause(0.001); end 示例 5：普通三维坐标下一个点沿着曲线移动的轨迹。 解： %三维螺旋线坐标 t1=10*pi*(0:;32 x1=cos(t1); y1=sin(t1); z1=-t1; %绘制曲线下面的水平直线 t2=(0:10)/10; x2=x1(end)*(1-t2); y2=y1(end)*(1-t2); z2=z1(end)*ones(size(x2)); %绘制垂直直线 t3=t2; z3=(1-t3)*z1(end); x3=zeros(size(z3)); y3=x3; %绘制曲线上面的水平直线 t4=t2; x4=t4; y4=zeros(size(x4)); z4=y4; %绘制曲线 x=[x1,x2,x3,x4]; y=[y1,y2,y3,y4]; z=[z1,z2,z3,z4]; plot3(x,y,z,'b','linewidth',3) axis off %绘制移动的点 h=line('Color',[1 0 0],'Marker','.','MarkerSize',40,'EraseMode','xor'); n=length(x); i=1;j=1; while 1 set(h,'xdata',x(i),'ydata',y(i),'zdata',z(i)); drawnow %使 Matlab 暂停目前的任务序列而去刷新屏幕 pause(0.0005); f=getframe(gcf); i=i+1; j=j+1; if i&n & j&450 end end33 5、Matlab 在高等数学中的应用5.0 Matlab 中的符号表示1、建立符号变量 Matlab 提供了 sym 和 syms 两个建立符号对象的函数 1）sym 函数，用来建立单个符号变量，调用格式为： 符号变量名=sym(?符号字符串?) 注：符号字符串可以是 常量、变量、函数或者表达式。 2）syms 函数，依次可以定义多个符号变量，调用格式为： syms 符号变量 1 符号变量 2 …… 符号变量 n 注意 符号变量名之间不能使用任何标点符号，只能用空格隔开。 例如： f=sym(?cos(x)?),f=sym(?sin(x)^2=0?); f=sin(x)+cos(x), 2、建立符号表达式 ? ? ? 利用单引号来生成表达式 利用 sym 函数建立符号表达式 使用已经定义过的符号变量组成符号表达式 注意：符号表达式包括符号函数和符号方程，区别在于是否带有等号。 3、符号表达式的运算 (1)四则运算 四则运算分别使用 symadd, symsub,symmul 和 symdiv 来实现加减乘除运算，使用 sympow 实现幂运算。 （注：6.5 以后的版本中已经没有了这些函数，而直接使用+ - * / 即可） (2)提取分子和分母 如果符号表达式是一个有理分式或可以展开为有理分式， 可利用 numden 函数来提 取符号表达式中的分子或者分母，调用格式为：34 [n,d]=numden(s) [n,d] = numden(sym(4/5)) returns n = 4 and d = 5. [n,d] = numden(x/y + y/x) returns (3)因式分解与展开 1）factors(s)分解因式 2）expand(s)展开 3）collect(s)合并同类项 4）collect(s,v)按变量 v 进行合并同类项 (4)符号表达式的化简 1）simplify(S)应用函数规则对 S 进行化简 2）simple(S) 进行综合化简 (5) 符号表达式与数值表达式之间的转换 1）利用函数 sym 可以将数值表达式表示成符号表达式 2）numeric 或者 eval 函数可以将符号表达式转换成数值表达式 3）函数 digits 和 vpa 配合替换函数 subs 进行转换。 digits 函数，digits(D)函数设置有效数字个数为 D 的近似解精度。 vpa 函数 ，R=vpa(S)符号表达式 S 在 digits 函数设置下的精度的数值解。vpa(S,D) 符号表达式 S 在 digits(D)精度下的数值解。 subs 函数，subs(S, NEW, OLD) 4、符号函数的相关运算 （1）复合函数运算。 compose 函数 （2）反函数运算。 finverse 函数 5、符号代数方程求解 (1)线性方程组的求解，函数 linsolve，solve，可以得到方程的精确解 (2)多变量的非线性方程的符号解法，使用函数 fsolve，调用格式有： ? ? ? ? ? X=fsolve(?fun?,X0) X=fsolve(?fun?,X0,options)options 为选择参数输入向量 X=fsolve(?fun?,X0,options,?gradfun?),gradfun 为输入函数在 X 处的偏导数 X=fsolve(?fun?,X0,options,?gradfun?,P1,P2,…)P1,P2 为问题定性参数 [X,options]=fsolve(?fun?,X0,…)返回使用的优化方法的参数 n = x^2+y^2 , d = y*x35 注意：fun 必须使用@指定 示例：定义一个函数 function F = myfun(x,c) F = [ 2*x(1) - x(2) - exp(c*x(1)) -x(1) + 2*x(2) - exp(c*x(2))]; end 求解该函数： c = -1; %定义输入参数 x = fsolve(@(x) myfun(x,c),[-5;-5]) 6、多项式的表示方法 （1）多项式的表示方法――转化为向量问题 对于多项式P( x) ? a0 x n ? a1 x n ?1 ? ... ? an?1 x ? an用行向量表示： P ? [a , a ,..., a , a ] 0 1 n ?1 n ①系数向量直接输入法， MATLAB 自动将向量元素按降幂顺序分配给各系数值。 函数 poly2sym 可以将向量表示的多项式转化为符号多项式表示。 注：由特征多项式生成 的多项式首项系数一定为 1；n 阶矩阵一般产生 n 次多项式。 ③由根创建多项式，由函数 poly 实现。注：若要生成实系数多项式，则根中的复数必 定对应共轭；生成的多项式向量包含很小的虚部时可用 real 命令将其过滤掉。 （2）多项式的运算 ①多项式求值。 输入变量值代入多项式计算时以数组为单元的使用函数 polyval（对应 元素计算）;以矩阵(必须为方阵)为计算单元求多项式的值用函数 ②多项式求根。两种方法，一种是调用函数 roots，另一种是通过建立多项式的伴随矩 阵再求其特征值的方法得到多项式的所有根。 （使用 compan 和 eig 函数） ③多项式的乘除法运算。乘法使用函数 conv（向量卷积） ，除法使用函数 deconv ④多项式微分。微分函数 polyder ⑤多项式拟合。两种方法，一种是由矩阵的除法求解超定方程来进行，另一种是用拟 合函数 polyfit,调用方式为 polyfit(X,Y,n)和[p,s]=polyfit(X,Y,n)36 5.1 求极限Matlab 求极限的命令为 limit(expr, x, a) limit(expr, a) limit(expr) limit(expr, x, a, 'left') limit(expr, x, a, 'right') 其中 limit(expr, x, a)表示求符号表达式 expr 关于符号变量 x 趋近于 a 时的极限， limit(expr)求表示缺省变量趋近于 0 时的极限。 例 15 求下列表达式的极限解 （1）clc, clear syms x f=(3*sin(x)+x^2*cos(1/x))/((1+cos(x))*log(1+x)); s=limit(f,x,0) s=simplify(s)（2）clc, clear syms x a s1=limit(((x+2*a)/(x-a))^x,x,inf) s2=limit(((x+2*a)/(x-a))^x,x,-inf)5.2 求导数Matlab 的求导数命令为 diff(expr) diff(expr, v) diff(expr, sym('v')) diff(expr, n) diff(expr, v, n)37 diff(expr, n, v) 其中 diff(expr)表示求表达式 expr 关于默认变量的 1 阶导数，diff(expr, v, n)和 diff(expr, n, v)都表示求表达式，expr 关于符号变量 v 的 n 阶导数。解 （1）clc, clear syms x y=log(x+sqrt(1+x^2)); s=diff(y,x,2) s=simple(s) %对符号函数进行化简 （2）a=[1,0.5,3.5,6]; da=diff(a)5.3 求一元函数极值计算的 Matlab 程序如下 syms x y=x^3+6*x^2+8*x-1; dy=diff(y); dy_zero=solve(dy), dy_zero_num=double(dy_zero) %变成数值类型 ezplot(y) %符号函数画图5.4 求积分1. 求不定积分 Matlab 求符号函数不定积分的命令为 int(expr) int(expr, v)38 解 syms x I=int(1/(1+sqrt(1-x^2))); pretty(I) %写成符号表达式的形式 2. 求定积分 （1）求定积分的符号解 Matlab 求符号函数的定积分命令为 int(expr, a, b) int(expr, v, a, b)解 syms x I=int(cos(x)*cos(2*x),-pi/2,pi/2) （2）求定积分的数值解 例 20 求下列积分的数值解输入 I=quadl(@(t) (t-3*t.^2+2*t.^3).^(-1/3),eps,0.5) 得到 I=1.4396。 注：quad()函数使用 simpson 法，即用二次曲线逼近被积函数，quadl()采用更高阶 的逼近方法。它们的调用方法是：q=quad(fun,a,b,tol)或者 quadl(fun,a,b,tol)，fun 可以是39 符号方程，也可以是匿名函数。 ）输入 I=dblquad(@(x,y)sqrt(1-x.^2-y.^2).*(x.^2+y.^2&=x),0,1,-0.5,0.5) 得到 I= 0.6028。 注： Matlab 中计算二重积分的函数形式为 q=dblquad(fun,xmin,xmax,ymin,ymax,tol)， 一般只用于积分区域为矩形的情况。为了使用该函数，有时需要将非矩形区域化成矩形 区域。2 2 示例：计算二重积分 ?? ?( x ? y )dxdy 积分区域 ? 由 x=1,y=x 及 y=0 围成。 ??解：首先绘制积分区域，Matlab 程序如下： clear,clc fill([0,1,1,0],[0,0,1,0],'y');hold on%绘制积分区域 fill([0.55,0.6,0.6,0.55,0.55],[0,0,0.6,0.55,0],'r') %绘制单元条 gtext('y=x'); gtext('x=1');pause gtext('y=0') 按照矩形区域调用 dblquad 函数，程序如下： I=dblquad(?(x.^2+y.^2).*(y-x&0)?,0,1,0,1) 求得 I=0.3333 （3）输入 fun3=@(x,y,z)z.^2*log(x.^2+y.^2+z.^2+1)./(x.^2+y.^2+z.^2+1).*(z&=0 &... z&=sqrt(1-x.^2-y.^2)); %该语句使用了续行符... I=triplequad(fun3,-1,1,-1,1,0,1) 得到 I= 0.1273。 %三重积分计算函数40 Matlab 符号求解的程序如下 clc, clear syms r z theta I1=int((r^2+z^2)*r, z, r^2,sqrt(2-r^2)); %求最内层积分 I2=int(I1,r,0,1); %求中间层的积分 I=int(I2,theta,0,2*pi) %求最终的积分结果 pretty(I) %分数线中间显示的格式 Inum=double(I) %把符号解化成数值解 解法二 求数值解的 Matlab 程序如下 clc, clear f=@(x,y,z) (x+y+z).^2.*(z&=x.^2+y.^2 & x.^2+y.^2+z.^2&=2); I=triplequad(f,-sqrt(2),sqrt(2),-sqrt(2),sqrt(2),0,sqrt(2)) 求得积分的值为 2.4486。 注：quad 与 int 的区别，int 可以求具有解析解的不定积分和定积分，quad 只能求解 定积分。5.5 求解微分方程（组）1、微分方程(组)的解析解41 在 MATLAB 中用 D 表示导数， 例如 Dy 表示 y?， D2y 表示 y??， Dy(0)=5 表示 y?(0)=5 等，符号常微分方程求解通过函数 dsolve 实现，调用格式为： dsolve(?e?,?c?,?v?)包括微分方程，初始条件和指定变量三个部分，求解常微分方程 e 在初始条件 c 下的特解，v 描述方程中的自变量，省略则默认以 t 为自变量，若没有初 始条件，则获得通解。 例 1：求 du ? 1 ? u 2 的通解 dt解：Matlab 程序为： dsolve('Du=1+u^2','t') ?d2 y dy ? 4 ? 29 y ? 0 ? 的特解 例 2：求微分方程 ??dx 2 dx ? ??y (0) ? 0, y?(0) ? 15 解：Matlab 程序为： dsolve('D2y+4*Dy+29*y=0','y(0)=0,Dy(0)=15','x') ? dx ? dt ? 2 x ? 3 y ? 3 z ?? ? dy 例 3：求常微分方程组的通解 ? ? 4 x ? 5 y ? 3 z ??dt ? dz ? ? 4x ? 4 y ? 2z ??dt 解：Matlab 程序为： [x,y,z]=dsolve('Dx=2*x-3*y+3*z','Dy=4*x-5*y+3*z','Dz=4*x-4*y+2*z') 一阶齐次线性常微分方程组和非齐次线性方程组可以直接通过矩阵操作求解。 对于一阶齐次线性微分方程组： X ? ? AX , X ? ? x1 ,xn ? , A ? (aij )n?n 其解的形式为：X (t ) ? e A(t ?t0 ) X 0对于非齐次线性方程组： X ? ? AX? f (t ), X (t 0 ) ? X 0 解的形式为：X (t ) ? e A(t ?t0 ) X 0 ???? e A(t ? s ) f (s)ds? ? ? ?? 例 4： X ? ? ? 0 2 ?1? X , X (0) ? ?2 ?? ? ?0 0 3 ? ?1 ?2 1 3? ?1 ??42 解：a=[2,1,3;0,2,-1;0,0,2]; x0=[1;2;1];x=expm(a*t)*x0 pretty(x); %化简结果?1 0 0 ? ? 0 ? ? 0 ?? ? ? ? ? ? ?? 例 5： X ? ? ? 2 1 ?2 ? X ? ? 0 ? , X (0) ? ?1 ?? t ? ? ? ?3 2 1 ?e cos(t ) ?1解： a=[1 0 0;2 1 -2;3 2 1]; fs=[0;0;exp(s)*cos(2*s)]; x0=[0;1;1]; tx=int(expm(a*(t-s))*fs,s); %先求不定积分，再计算定积分的结果，提高速度 xstar=subs(tx,s,t)-subs(tx,s,0); x=expm(a*t)*x0+ x=simple(s),pretty(x) %simple 直接显示化简结果，simplify 显示化简过程，pretty 修 改显示格式 2、微分方程（组）的数值解 当难以求得微分方程的解析解时，可以求其数值 解。Matlab 中求解微分方程数值解 的方法有以下几个：ode45，ode23，ode113，ode23s，ode15s。调用形式如下： [t,x]=solver(‘f’, ts, x0, options)其中，t 是自变量的值，x 是函数值，solver 表示上述 5 种求解器的一种；‘f’是由待解方程写成的 m 文件名，ts=[t0, tf]表示自变量的初值和 终值，options 用于设定误差限，通过函数 options=odeset(‘reltol’, rt, ‘abstol’, at)分别设 定相对误差和绝对误差。 不同的函数代表不同的内部算法，ode45 表示 4/5 阶龙格-库塔-费尔贝格算法（RK 算法），是解非刚性常微分方程的首选方法，ode23 采用 2/3 阶 RK 方法，ode113 采用 多步法，效率一般比 ode45 高。ode15s，ode23s,ode23t,ode23tb 用来解刚性常微分方程。 注 1：在解 n 个未知函数的的方程组时，x0 和 xn 均为 n 为向量，m 文件中的待解 方程组应以 x 分量形式写出； 注 2：使用 Matlab 求数值解时，高阶微分方程必须等价的变换成一阶微分方程组。 ? ?d2 x 2 dx ? ? ? 1000(1 ??x ) x 0 ? 例 1：求解 van der Pol 刚性微分方程 ??dt 2 的数值解 dt ? ??x(0) ? 2; x '(0) ? 043 解：令 y1=x，y2=y1?，则微分方程变为一阶微分方程组： ? y1' ? y 2 ? 2 ??y 2 ' ? 1000(1 ? y1 ) y 2 ? y1 ? y1(0) ? 2; y 2(0) ? 0 ?? 写成 m 文件如下所示： function dy=fun1(t,y) dy=zeros(2,1); dy(1)=y(2); dy(2)=1000*(1-y(1)^2)*y(2)-y(1); 取 [t0 ,tf]=[0,3000] ，输入命令： [t,y]=ode15s('fun1',[0,3000],[2,0]); plot(t,y(:,1),'-') ? y1 ' ? y2 y3 ?? ? ? y1 y3 ? y2 ' 例 2：求解微分方程组 ?? ? y3 ' ? ?0.51 y1 y2 y1 (0) ? 0, y2 (0) ? 1; y3 (0) ? 1 解：建立方程组的 m 文件： function dy=fun2(t,y) dy=zeros(3,1); %初始化 dy(1)=y(2)*y(3); dy(2)=y(1)*y(3); dy(3)=-0.51*y(1)*y(2); 取 t0=0，tf=12，Matlab 命令如下： [t,y]=ode45('fun2',[0 12],[0 1 1]) plot(t,y(:,1),'-',t,y(:,2),'*',t,y(:,3),'+')?? x ? ? ( y ??x) ?? y ? ? x ? y ? xz 在初始条件[x0,y0,z0]=[5,13,17]时的解，其 例 3：求解 Lorez 方程 ?? ?? z ? ?? z ??xy ??中 ? ? 10, ? ? 28, ? ? 8 / 3 。 解：rho=10;beta=29;lamda=8/3;44 f=@(t,Y)[rho*(Y(2)-Y(1));beta*Y(1)-Y(2)-Y(1)*Y(3);-lamda*Y(3)+Y(1)*Y(2)];%方程 组必须使用列向量表示 [t,y]=ode45(f,[0,30],[5,13,17]) subplot(2,2,1) plot(t,y(:,1),'*') %x 曲线 subplot(2,2,2) plot(t,y(:,2),'X') subplot(2,2,3) plot(t,y(:,3),'O') subplot(2,2,4) plot3(y(:,1),y(:,2),y(:,3)) %绘制空间轨迹图 3、边值问题的 Matlab 数值解 Matlab 中使用 bvp4c（bvp5c）和 bvpinit 命令求解常微分方程中的两点边值问题， 其调用格式为： sol = bvp4c(odefun,bcfun,solinit,options) solinit = bvpinit(x,yinit,parameters) %该函数用来给 sol 提供初始猜测解。5.6 级数求和Matlab 级数求和的命令为 r = symsum(expr, v) r = symsum(expr, v, a, b) 其中 expr 为级数的通项表达式，v 是求和变量，a 和 b 分别为求和变量的起始点和 终止点，若没有指明 a 和 b，a 的默认值为 0，b 的默认值为 v-1。 无穷数列的累加称为级数，当取其前面若干有限项时，得到的是部分和，数列 a 累 加形成的新序列可用 s=cumsum(a)实现。 a 的长度为 n， 则 s 的长度也为 n， 即每一个 s(k) 是数组 a 中前 k 项的和。注意 cumsum 与 sum 的区别，sum(a)得到的是一个数，即序列 s 中最后一项。 示例： clc,clear k=1:10; a=1./k.^2;45 s=cumsum(a) ss=sum(a) 例 22 求如下级数的和解 （1）syms n f1=(2*n-1)/2^n; s1=symsum(f1,n,1,inf) （2）syms n f2=1/n^2; s2=symsum(f2,n,1,inf) 补充一点： Matlab 求数组元素乘积的函数 B=prod(A) 如果 A 是一个 m 行一列的 （向 量），则这种用法即返回这 m 个元素的乘积； 如果 A 是一个 m 行 n 列的矩阵，则 A 的每一列都被看做一个 m 行 1 列的向量，分别计算每个向量中元素的乘积，返回给 B， 因此 B 是一个 1 行 n 列的数组。5.7 Matlab 求解优化问题Matlab 求解各种形式的极值问题被集成在优化工具箱中。 下面是优化工具箱 4.0 版本的功能示意图：46 47 变量 f描述 线性规划目标函数 f*X 或二次规划的目标 函数 X’*H*X+F*X 中线性项的系数向量 非线性优化的目标函数. fun 必须为行命令 对象或 M 文件、 嵌入函数、 或 MEX 文件的 名称 二次规划目标函数 X’*H*X+f*X 中的二次 项的系数矩阵 A 矩阵和 b 向量分别为线性不等式约束： AX≤b 中的系数矩阵和右端向量 Aeq 矩阵和 beq 向量分别为线性等式约束 Aeq*X=beq 中的系数矩阵和右端向量 X 的下限和上限向量 迭代初始点坐标 函数最小化的区间 优化选项参数结构调用函数 linprog,quadprog fminbnd,fminsearch,fminunc, fmincon,lsqcurvefit,lsqnonlin, fgoalattain,fminimax quadprog linprog,quadprog,fgoalattain, fmincon,fminimax linprog,quadprog,fgoalattain, fmincon,fminimax linprog,quadprog,fgoalattain, fmincon,fminimax,lsqcurvefit, lsqnonlin 除 fminbnd 外所有函数 fminbnd 所有优化函数funHA,bAeq,beqvlb,vub X0 x1,x2 options变量描述 描述退出条件： exitflag&0 表示目标函数收敛于 x exitflag=0 表示已达到函数评价或迭代的最大次 数 exitflag&0 表示目标函数不收敛 由优化函数求的 x 值调用函数exitflagx所有优化函数 linprog,quadprog,fgoalattain fmincon,fminimax,lsqcurvefit lsqnonlin,fminbndfeval解 x 处的目标函数值output包含优化结果信息的输出结构 所有优化函数 iterations,Algorithm,FuncCount（函数评价次数）进行极值计算时有两种方式，通过命令 optimtool 打开优化工具箱 GUI 求解和使用48 上述提供的命令函数进行求解。 示例一：使用单纯形计算最小值 min f ? ?4 x1 ? x2 ? ? x1 ? 2 x2 ? 4 ? 2 x ? 3 x ? 12 ?? 2 s.t. ? 1 ? x1 ? x2 ? 3 x1 , x2 ? 0 方法一：通过优化工具箱求解 方法二：输入 Matlab 命令： clc,clear f=[-4;-1]; A=[-1 2;2 3;1 -1]; b=[4;12;3]; [x,fval,exitflag,output]=linprog(f,A,b) 示例二：非线性优化 f ( x) ? ? x1 x2 x3 ? x ? 2 x2 ? 2 x3 ? 0 ?? s.t. ? 1 ? x1 ? 2 x2 ? 2 x3 ? 72 方法一：通过优化工具箱求解 方法二：输入 Matlab 命令： clc,clear f=@(x)(-x(1)*x(2)*x(3)); A=[-1 -2 -2;1 2 2]; b=[0;72]; [x,fval,exitflag,output]=fmincon(f,[10;10;10],A,b) 其他求解极值的函数主要 有： (1)x=fminbnd(fun,x1,x2)，求函数在区间[x1,x2]上的极小值对应的自 变量值 示例三：求 f(x)=x^4-x^2+x-1 在区间[-2,1]上的极小值。 解： Matlab 命令如下： [x,feval,exitflag,output]=fminbnd('x^4-x^2+x-1',-2,1) (2) x=fminsearch(fun,x0)无约束非线性极小值求解，x0 为起始搜索点。49 示例四：fminsearch 函数求 f(x)=1/((x1-2)^2+3)-1/(2(x2+1)^2-5),初始点[0,0]， (3) x=fminunc(fun,x0,options)求解函数的一个局部极小值 (4) x=fminimax(fun,x0,A,b,Aeq,beq,lb,ub,nonlcon)nonlcon 包括非线性约束; (5) x=lsqnonlin(fun,x0,lb,ub);非线性最小二乘优化 (6) x=linprog(f,A,b,Aeq,beq,lb,ub,x0)求解线性规划问题 (7) x=bintprog(f,A,b,Aeq,beq,x0)求解 0-1 规划 (8) x=quadprog(H,f,A,b,Aeq,beq,lb,ub,x0);求解二次规划5.8 Matlab 求解插值和拟合问题5.8.1 Matlab 插值工具箱 1、一维插值函数 interp1，命令格式为： yi=interp1(x,Y,xi,?method?)其中 method 指定的方法为 nearest,linear,spline,cubic， pchip （分段三次 Hermite 插值） ， 所有的插值方法都要求 x 是单调的。 x 与 Y 是具有相同大小 的向量，求在 xi 点处的插值函数值 yi 一维快速插值函数 interp1q，命令格式为： yi = interp1q(x,Y,xi)，其中 x 是单调递增的列向量，Y 是列向量或者行数为 length(x) 的矩阵，xi 是一个列向量。 yi = interp1(x,Y,xi,method,'extrap') %对于超出 x 范围的 xi 中的分量将执行特殊的外 插值法 extrap。 示例 1：在一 天 24 小时内，从零点开始每间隔 2 小时测得的环境温度数据分别为 12，9，9，10，18 ，24，28，27，25，20，18，15，13，推测中午 12 点（即 13 点）时的温 度． 解： x=0:2:24; y=[12 9 9 10 18 24 28 27 25 20 18 15 13]; a=13; y1=interp1(x,y,a,'spline') xi=0:1/3600:24; %绘制一天内的温度变化曲线 yi=interp1(x,y,xi, 'spline'); plot(x,y,'o' ,xi,yi) 示例 2： x=0:0.12:1; y=(x.^2-3*x+5).*exp(-5*x).*sin(x); plot(x,y,'ro',x,y);50 x1=0:0.02:1; %要插值点 y1=(x1.^2-3*x1+5).*exp(-5*x1).*sin(x1); y2=interp1(x,y,x1);y3=interp1(x,y,x1,'spline'); y4=interp1(x,y,x1,'nearest');y5=interp1(x,y,x1,'cubic'); plot(x1,[y2' y3' y4' y5'],':',x,y,'ro',x1,y1); legend('linear','spline','nearest','cubic','样本点','原函数'); %计算各插值的最大误差 [max(abs(y1-y2)),max(abs(y1-y3)),max(abs(y1-y4)),max(abs(y1-y5))]; 2、二维数据内插值 interp2，插值函数为二元函数。 ZI = interp2(X,Y,Z,XI,YI,?method?) X 和 Y 分别是 m 维和 n 维的向量，表示节点，Z 为 n*m 的矩阵，表示节点值，XI， YI 是为一维数组，表示插值点，XI 与 YI 是方向不同的向量，一个是行向量，一个是列 向量，ZI 是行数为 YI 的维数，列数为 XI 的维数的矩阵，表示得到的插值。若 XI 与 YI 中有在 X 与 Y 范围之外的点，则相应地返回 nan（Not a Number） 。method 为插值 方法。 如果在某区域测量了若干个节点的高程（节点值） ，为了画出比较精确的等高线图， 需要先插入更多的点（插值点） ，然后计算这些点的高程。 如果是三次样条插值，可以使用命令 pp=csape({x0,y0},z0,conds,valconds);z=fnval(pp,{x,y}) 示例 1： years = 0; service = 10:10:30; wage = [150.697 199.592 187.625 179.323 195.072 250.287 203.212 179.092 322.767 226.505 153.706 426.730 249.633 120.281 598.243]; w = interp2(service,years,wage,15,1975) %计算在 1975 年工作了 15 年的员工的工资[x,y]=meshgrid(-3:0.6:3,-2:0.4:2); z=(x.^2-2*x).*exp(-x.^2-y.^2-x.*y); surf(x,y,z);figure [x1,y1]=meshgrid(-3:0.2:3,-2:0.2:2); z1=interp2(x,y,z,x1,y1); z2=interp2(x,y,z,x1,y1,'cubic'); z3=interp2(x,y,z,x1,y1,'spline'); surf(x1,y1,z1);figure %比较误差 z0=(x1.^2-2*x1).*exp(-x1.^2-y1.^2-x1.*y1);51 surf(x1,y1,abs(z0-z2));axis([-3,3,-2,2,0,0.08])surf(x1,y1,abs(z0-z3)); axis([-3,3,-2,2,0,0.08]) interp2 函数能够较好的进行二维插值运算，但是它只能处理以网格形式给出的数 据，如果数据是随机给出的，此时可以使用 griddata 函数。其调用格式为： 格式:z=griddata(x0,y0,z0,x,y,?method?) ? v4 ?:MATLAB4.0 提供的插值算法，公认效果较好； ?linear ?：双线性插值算法（缺省算法） ； ?nearest?：最临近插值； ?cubic?：双三次插值。clear x=-3+6*rand(199,1);y=-2+4*rand(199,1); z=(x.^2-2*x).*exp(-x.^2-y.^2-x.*y); plot(x,y,'*');plot3(x,y,z,'*'); [x1,y1]=meshgrid(-3:0.2:3,-2:0.2:2); z1=griddata(x,y,z,x1,y1,'cubic'); surf(x1,y1,z1); z2=griddata(x,y,z,x1,y1,'v4');surf(x1,y1,z2); %误差比较 z0=(x1.^2-2*x1).*exp(-x1.^2-y1.^2-x1.*y1); surf(x1,y1,abs(z0-z1));axis([-3,3,-2,2,0,0.15]) surf(x1,y1,abs(z0-z2));axis([-3,3,-2,2,0,0.15]) 3、 三维插值运算函数 interp3,n 维网格插值 interpn， 其调用格式同 interp1 和 interp2， 对应的三维网格生成函数为[x,y,z]=meshgrid(x1}
&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&签到排名：今日本吧第个签到，本吧因你更精彩，明天继续来努力！
本吧签到人数：0成为超级会员，使用一键签到本月漏签0次！成为超级会员，赠送8张补签卡连续签到：天&&累计签到：天超级会员单次开通12个月以上，赠送连续签到卡3张
关注：172,543贴子：
关于matlab解微分方程的问题，函数含有未知常数，怎么定义函数？收藏
请问各位大神，如果我想把这个方程F(q) = dq/dt = q(1 - q)( pr4 - c2) F( p) = dp/dt = p(1 - p)(qr1 + qc3 - qr2 - c1)在matlab里面定义出来，然后再画出它的相位图，其中p,q是未知变量，其他的r1,r2,r4,c1,c2,c3是未知常数，最后求出的解中可以用其作为表达式，我该怎么写函数，写的函数好像都有问题，有大神帮忙吗？我写了两种，好像都通不过，大神帮忙看一下：一种：function dxdt=eq2(t,x)global a b c d e fdxdt=[x(1)*(1-x(1))*((a-b+e)*x(2)-d);x(2)*(1-x(2))*(c*x(1)-f)];end另一种：function xdot=eq3(t,x)%创建数组储存数据global a b c d e fxdot=zeros(2,1);xdot(1)=x(1)*(1-x(1))*((a-b+e)*x(2)-d);xdot(2)=x(2)*(1-x(2))*(c*x(1)-f);
网站建设10年沉淀，搭建专属品牌网站
楼主我也遇到这个问题了，微分方程里面有未知常数，这种方程怎么求解
登录百度帐号推荐应用查看: 30172|回复: 18|关注: 1
matlab dsolve求解微分方程组
<h1 style="color:# 麦片财富积分
新手, 积分 25, 距离下一级还需 25 积分
关注者： 5
最近有人问有关微分方程求解问题 特举一例供大家参考：
dsolve调用格式：dsolve（'equ1','equ2',..........'equN'）
另外要注意：在微分方程表达式输入中，以大写字母D来表示微分，D,D2,.......Dn分别表示一阶，二阶和n阶
2dx/dt+dy/dt-y=exp(-t)
dx/dt+x+y=0& && && && && && && && && &&&其中初始条件：x(0)=1.5,y(0)=0
首先求解微分方程的通解：
s=dsolve('2*Dx+Dy-y=exp(-t)','Dx+x+y=0');%求解微分方程组的通解
%微分方程组变量x的通解
-C1*exp((1+2^(1/2))*t)-C2*exp(-(2^(1/2)-1)*t)+1/2*C1*exp((1+2^(1/2))*t)*2^(1/2)-1/2*C2*exp(-(2^(1/2)-1)*t)*2^(1/2)-1/2*exp(-t)
%微分方程组变量y的通解
C1*exp((1+2^(1/2))*t)+C2*exp(-(2^(1/2)-1)*t)
然后根据初始条件，求解微分方程组的特解：
&& s=dsolve('2*Dx+Dy-y=exp(-t)','dx+x+y=0','x(0)=1.5','y(0)=0');%微分方程组在给定初始条件下的特解
Warning: Explicit solution could not be found.
& In dsolve at 333
&& s=dsolve('2*Dx+Dy-y=exp(-t)','Dx+x+y=0','x(0)=1.5','y(0)=0');%微分方程组在给定初始条件下的特解
-2^(1/2)*exp((1+2^(1/2))*t)+2^(1/2)*exp(-(2^(1/2)-1)*t)+exp((1+2^(1/2))*t)+exp(-(2^(1/2)-1)*t)-1/2*exp(-t)
2^(1/2)*exp((1+2^(1/2))*t)-2^(1/2)*exp(-(2^(1/2)-1)*t)
&& %或者使用下面命令直接获取x,y的特解
&& [x,y]=dsolve('2*Dx+Dy-y=exp(-t)','Dx+x+y=0','x(0)=1.5','y(0)=0')
-2^(1/2)*exp((1+2^(1/2))*t)+2^(1/2)*exp(-(2^(1/2)-1)*t)+exp((1+2^(1/2))*t)+exp(-(2^(1/2)-1)*t)-1/2*exp(-t)
2^(1/2)*exp((1+2^(1/2))*t)-2^(1/2)*exp(-(2^(1/2)-1)*t)
[ 本帖最后由 edifiers2008 于
11:12 编辑 ]
<h1 style="color:# 麦片财富积分
论坛优秀回答者
帖子最佳答案
关注者： 595
如果用simulink怎么求解非线性微分方程组呢?
<h1 style="color:# 麦片财富积分
<h1 style="color:# 麦片财富积分
请教您一个问题，为什么我的MATLAB RA2009中的dsolve函数不能调用呢？使用help dsolve命令说没有找到，用dsolve求解时就会报错说没有定义的函数。谢谢
<h1 style="color:# 麦片财富积分
我的、dsolve命令也不可用也是2009a
<h1 style="color:# 麦片财富积分
自己在网上找一个，拷进去就可以用
<h1 style="color:# 麦片财富积分
谢谢楼主！学习了，哈哈！
<h1 style="color:# 麦片财富积分
回复 6# damaomao 的帖子
2009a的dsolve命令不可用是怎么解决的？谢谢
<h1 style="color:# 麦片财富积分
新手，学习一下！
站长推荐 /2
快速搭建新能源汽车整车模型及其性能优化
MATLAB中文论坛是全球最大的 MATLAB & Simulink 中文社区。用户免费注册会员后，即可下载代码，讨论问题，请教资深用户及结识书籍作者。立即注册加入我们吧！
MATLAB官方社交平台
MATLAB中文论坛微社区4084人阅读
Matlab（10）
自动控制原理（3）
&[time,x]=solver(str,t,x0)&& 计算ODE或由字符串str 给定的ODE的值，部分解已在向量time中给出。在向量time
中给出部分解，包含的是时间值。还有部分解在矩阵x中给出，x的列向量每
个方程在这些值下的解。对于标量问题，方程的解将在向量x中给出。这些解
在时间区间t(1)到t(2)上计算得到。其初始值是x0 即x(t(1)).此方程组有str 指定的
M文件中函数表示出。这个函数需要两个参数：标量t 和向量x, 应该返回向量
x’(即x的导数)。因为对标量ODE来说，x和x’都是标量。在M文件中输入odefile
可得到更多信息。同时可以用命令numjac来计算Jacobi函数。
[t,x]=solver(str,t,x0,val)&& 此方程的求解过程同上，结构val包含用户给solver的命令。参见odeset和表1，可得到更多信息。
Ode45&&&&&此方法被推荐为首选方法。
Ode23&&&&&这是一个比ode45低阶的方法。
Ode113&&&&&用于更高阶或大的标量计算。
Ode23t&&&&&用于解决难度适中的问题。
Ode23s&&&&&用于解决难度较大的微分方程组。对于系统中存在常量矩阵的情况也有用。
Ode15s&&&&&与ode23相同，但要求的精度更高。
Ode23tb&&&&&用于解决难度较大的问题，对于系统中存在常量矩阵的情况也有用。
Set=odeset(set1,vak1,set2,val2,…)& 返回结构set,其中包含用于ODE求解方程的设置参数，&& 有关可用设置
的信息参见表1。
Odeget(set,’set1’)&&& 返回结构set中设置set1的值。
有许多设置对odeset控制的ODE解是有用的，参见表1。例如，如果要在求解过程中画出解的图形，可以
输入：inst=odeset(‘outputfcn’,’odeplot’);.也可使用命令odedemo。
表1&&& ODE求解方程的设置参数
RelTol&&&给出求解方程允许的相对误差
AbsTol&&给出求解方程允许的绝对误差
Refine&&给出与输入点数相乘的因子
OutputFcn&这是一个带有输入函数名的字符串，该字符串将在求解函数执行的每步被调用：odephas2(画出2D
的平面相位图)。Odephas3( 画出3D 的平面相位图)，odeplot( 画出解的图形)，odeprint( 显示中间结
OutputSel&是一个整型向量。指出哪些元素应该被传递给函数，特别是传递给OutputFcn
Stats&&如果参数Stats为on,则将统计并显示出计算过程中资源消耗情况
Jacobian… 如果编写ODE文件代码以便F（t,y,jocobian）返回dF/dy,则将jacobian设置为on
Jconstant…如果雅可比数df/dy是常量，则将此参数设置为on
Jpattern…& 如果编写ODE文件的编码以便函数F（[],[],jpattern）返回带有零的稀疏矩阵并输出非零元素dF/fy,则
需将Jpattern设置为on
Vectorized…如果编写ODE文件的编码以便函数F(t,[y1,y2……])返回[F(t,y1)& F(t,y2)…],则将此参数设置成on
Events…&& 如果ODE文件中带有参数‘events’，则将此参数设置成on
Mass…&& 如果编写ODE文件编码以实现函数F(t,[],‘mass’)返回M和M(t),应将此参数设置成on
MassConstant…如果矩阵M(t)是常量，则将此参数设置成on
MaxStep…&& 此参数是限定算法能使用的区间长度上限的标量
InitialStep…&& 给出初始步长的标量。如果给定的区间太大，算法就使用一个较小的步长
MaxOrder…&& 此参数只能被ode15s使用，它主要是指定ode15s的最高阶数,并且此参数应是从1到5的整数
BDF…&& 此参数只能被ode15s使用，如果倒推微分公式而不是用通常所使用的微分公式，则要将它设置为
NormControl…如果算法根据norm(e)&=max(Reltol*norm(y),Abstol)来步积分过程中的错误，则要将它设置成on
求解x’=-x^2 初值x(0)=1
创建函数xprim1,将此函数保存在M文件xprim1.m中：
function xprim=xprim1(t,x)
xprim=-x.^2;
然后调用MATLAB的ODE算法求解方程。然后画出解的图形：
[t,x]=ode45(‘xprim1’,[0,1],1);
plot(t,x,’-‘,t,x,’o’);
xlabel(‘time t0=0,tt=1’);
ylabel(‘x values x(0)=1’);
首先创建xprim2,将此函数保存在M文件xprim2.m中：
function xprim=xprim2(t,x)
xprim=x.^2;
然后调用MATLAB的ODE算法求解方程。然后画出解的图形：
[t,x]=ode45(‘xprim2’,[0,0.95],1);
注意：在MATLAB中计算出的点在微分绝对值大的区域内更密集些。
如果想在图像中直接观察，则在代码后输入以下命令
plot(t,x,’o‘,t,x,’-’);
xlabel(‘time t0=0,tt=0.95’);
ylabel(‘x values x(0)=1’);
例2求解下列二元一阶微分方程组
X1’=X1-0.1X1X2+0.01t
X2’=-X2+0.02X1X2+0.044
这个方程组用在人口动力学中。可以认为是单一化的捕食者---被捕食者模式。例如，狐狸和兔子。表示被捕食者， 表示捕食者。如果被捕食者有无限的食物，并且不会出现捕食者。于是有，这个式子是以指数形式增长的。大量的被捕食者将会使捕食者的数量增长；同样，越来越少的捕食者会使被捕食者的数量增长。而且，人口数量也会增长。
创建xprim3,将此函数保存在M文件xprim3.m 中：
function xprim=xprim3(t,x)
xprim=[x(1)-0.1*x(1)*x(2)+0.01*t; …
-x(2)+0.02*x(1)*x(2)+0.04*t];
然后调用一个ODE算法和画出解的图形：
[t,x]=ode45(‘xprim3’,[0,20],[30;20]);
plot(t,x);
xlabel(‘time t0=0,tt=20’);
ylabel(‘x values x1(0)=30,x2(0)=20’);
在MATLAB中,也可以根据x2函数绘制出x1图形,命令plot(x(:2),x(:1)) 可绘制出平面相位图.
例3 带参数的微分方程
X1’=a-(b+1)X1+X1^2X2
X2’=bX1+X1^2X2
方程由下面的M文件stiff1.m定义:
functionstiff=stiff1(t,x)
&& %变量不能放入参数表中
stiff=[0,0];&&& %stiff必须是一个冒号变量
stiff(1)=a-(b+1)*x1+x(1)^2*x(2);
stiff(2)=b*x(1)-x(1)^2*x(2);
下面的M文件给出一个比较困难的问题：
[t,X]=ode23(‘stiff1’，[0 10],[1 3]);
运行后得到的结果如下：
elapsed_time=72.1647
用专门解决复杂问题的解法ode23s,将得到较好的结果：
elapsed_time=1.0098
对于边界值问题，除了微分方程，还有边界处的值。在一维下这意味着至少有两个条件。
例4假设要研究一根杆的温度分布情况。这根杆一端的温度是T0,另一端的温度是T1。令y(x)表示这根杆的温度，函数f(x)表示加热源。从时间t=0开始，在相当长的时间内加热这根杆，直至达到平衡状态。这就是所谓的定常值或稳定状态。这个定常值可由下面的方程模型表示：
-Y’(X)=f(X),0&X&1
假设这根杆两端为：x=0和x=1。
假设在其两端又一根固定的柱子（或者可以看成是一个连接两个岛屿的桥），如图7所示。
令y(x)表示加载函数g(x)后弯曲的柱子。此问题需要有两个关于此柱子两端的边界条件。假设这根柱子非常牢
固的固定在墙上，即y在墙上的导数是0。可以得到下面的ODE,其中介绍了自然协调系统：
Y’’’(X)=g(x),0&X&1
Y(0)= Y’(0)= Y(1)= Y’(1)=0
解y(x)的导数可由有限的差分代替
如果用这些差分方程来代替ODE中的导数，就能得到一个所有未知的yi的方程组。其系数矩阵是一个有序区
间，此区间的宽度决定于这个微分方程的导数个数.
根据前面的温度模型的方程研究一下杆的温度分布，将所有的导数换成不同的差分并得到：
其中fj=f(xj)。为了简单起见，设M=6,即给定的y0和y6，而y1,y2,….y5 为未知变量。于是就有
注意，y0=T0和yM=T1必须移到方程组的右边。此时得到的矩阵是一个对角矩阵，其对角线上的元素为2，并且上一对角线和下一对角线上的元素为1。
下面解此问题的文件temperature.m。用户必须先给出分段数及f(x)(用点符号)，最后给出T0和T1。有关稀疏矩阵的更多信息参见其他资料。
%杆上的温度分布，用T0和T1分别表示两端温度
%这根杆放在x坐标的0和1 区间上，并被分成M个子区间，每个子区间的长度为1/M
%创建稀疏矩阵方程Ax=b并求解
%矩阵A对角阵，并以稀疏矩阵的形式存储
M=input(‘Give the number ofsubintervals (M):’);
Deltax=1/M;
xx=0:deltax:1;
funcStr= input(‘give f(x),the extra heatsource(e.g.,x.^3):’,’s’);
T0=input(‘Give y(0) (left): ’);
T1=input(‘Give y(1) (right): ’);
%构造对角阵和方程右边b
vectorOnes=ones(M-1,1);
A=spdiags([-vectorOnes,2*vectorOnes,-vectorOnes],[-1 0 1],M-1,M-1);
x=xx(2:end-1);
%x为区域内的值。
f=eval(funcStr);
%响应的f(x)的值。
b=deltax^2f;
b(1)=b(1)+T0; %对边界值x=0,x=1 进行特殊处理。
b(end)=b(end)+T1;
%解线性方程
%y在区间内：j=1:M-1.
y=[T0;y;T1];
%y在整个区间内：0&=x&=1.
%上面图形表示外部热源。
%下面图形表示杆上的热分布。
subplot(2,1,1);
plot(x,f);
title(‘External heatsource f(x).’,’FontSize’,14);
subplot(2,1,2);
plot(xx,y,’r’);
title(‘Tempearture distribution in a rod.’,’FontSize’,14);
%将区间分成等份，根据方程f(x)=在图中可以得到解。
例5 Matlab 中二阶常微分方程的数值解法
y''-3y'+2y=1;& y(0)=1;&y'(0)=0; 用数值解法求y(0.5)
odefun=@(t,x)[x(2);3*x(2)-2*x(1)+1];
[t,y]=ode45(odefun,[0:0.01:2],[1 0]);
y= dsolve('D2y-3*Dy+2*y=1','Dy(0)=0','y(0)=1');
y =exp(t) - exp(2*t)/2 + 1/2
&& feval(@(t)exp(t) - exp(2*t)/2 +1/2,0.5)
ans = 0.7896
例6求解范德普方程
首先是函数定义：
function xdot = myFcn (t, x)
xdot = [x(2) ; u(1-x(1)^2)*x(2)-x(1)]
其次是主程序：
% Solving options predetermined
options=odeset('OutputFcn',@odephas2,'reltol',1e-5);
% Solving vdp equation withdifferent ini conditions
for a = -3 : 0.1 : 3
b1=sqrt(25-a^2);
b2=-sqrt(25-a^2);
% Solve the problem
[t,y]=ode45(@myFcn, [0 20],[a b1], options);
[t,y]=ode45(@myFcn, [0 20],[a b2], options);
例7x1.^2-10x1+x2.^2+b=0
&&&& x1*x2.^2+x1-10x2+8=0
用MATLAB实现，编写非线性方程组的M文件fx1.m如下所示：
&& function y=fx1(x)
&&y(1)=x(1)*x(1)-10*x(1)+x(2)*x(2)+8;
&&y(2)=x(1)*x(2)*x(2)+x(1)-10*x(2)+8;
编写非线性方程组导数的M文件dfx1.m
& function y=dfx1(x)
&&&& y(1)=2*x(1)-10;
&&&& y(2)=2*x(2);
&&&& y(3)=x(2)*x(2)+1;
&&&& y(4)=2*x(1)*x(2)-10;
附一篇文章《用MATLAB求解非线性微分方程》
原始地址：http://blog.csdn.net/anhuixue/article/details/7560558
总结一下MATLAB中求解微分方程的思路和步骤。MATLAB中的帮助文件，内容实在是特别给力，简明的不能再简明了。可惜了是英文写的，即使是天天要面对英语的我，也颇为头疼。那么以下我将结合我对MATLAB中帮助文档的理解，提出我对求解最简单微分方程的方法的总结。以求解一个经典的非线性方程为例。
一求解一下这个方程，判断她是否稳定。要是稳定，那么她是否存在极限环;
二.一般的计算机求解方程的方法：
首先把该方程改写成一个规范的形式，一般使用状态空间表示法；
调用已有的算法进行求解；
最后对得出的结果进行处理，比如画图之类的。
三.输入待求解的方程。
MATLAB中关于输入输入待解方程的语句特别简单。需要先定义一个普通函数，函数名任意，姑且叫做myFcn，格式如下
function xdot = myFcn (t, x)&
需要注意的是，函数必须含有t, x两个参数，名称可以自己任意定。xdot是这个函数本身的返回值，只出现在这个函数内部，因此也可以任意定。
定义中的x被视为是一个列向量，x(i)表示列向量中的第i个分量。那么F函数的每一个分量即简单地用表达时给出即可。其中的自变量可以引用x(i)。以范德普方程为例：
xdot = [x(2) ; u(1-x(1)^2)*x(2)-x(1)]&
于是，这两句话便构成了待解函数。
四.调用MATLAB函数进行求解
通常人工求解微分方程需要知道初始值，计算机求解也不例外。另外，由于非线性方程一般只有数值解，故计算精度也可以调整。这些都是可以自己调整的参数。
调用MATLAB计算求解常微分方程的模式很简单，格式为：
[t, x] = ode45(@myFcn, [开始时间 结束时间], [初始值列向量],options)
注意到求解出来的结果是[t, x]即一堆数，所以需要我们进行后处理比如画图之类的。
ode45表示了MATLAB中的一种内置计算微分方程的算法，最为常用，出于娱乐目的，就先用这个。
@表示待求解的方程，如myFcn。
[开始时间 结束时间]表示求解的时间段，不必定义间隔。
[初始值列向量]就不用多说了。通常在求解之前定义一个变量x0列向量表示初始值，然后输入起来就方便多了。
options本身是一个变量。注意，她的名称不固定，但是他是一个以结构体为类型的变量。options很重要，稍后讨论。
五.进行后处理。在得到[t, x]后可以自己用plot神马的进行画图。
六.options的用法。
options在MATLAB帮助中是这样定义格式的：
&&&&&&options = odeset (“Name”, Value, “Name”, Value, …)
意思是，options这个变量要用到odeset这个函数来对他进行赋值。而odeset这个函数的参数必须是成对出现的：一个名称对应一个数值。
我们要用到的是这样的：
options= odeset(“OutputFcn”, @odephas2,“reltol”, 1e-5);
注意，odeset中的参数名称不可任意定，因为都是预定义好的。”OutputFcn”使用预定义的函数进行输出，而与定义的函数有好几种，使用@进行选择，我们选odephas2即画相平面图(phase portrait)。之后的那个参数表示相对误差的最大允许值，不说也明白。
有一个问题，用odephas2画出的图图不好看，因为上面打了好多圈圈。解决办法是找到该文件，修改其中所有的plot语句即可，把那个”-o”去掉就行了。
&&相关文章推荐
* 以上用户言论只代表其个人观点，不代表CSDN网站的观点或立场
访问：66045次
排名：千里之外
原创：28篇
转载：13篇
(3)(1)(1)(6)(9)(14)(7)
(window.slotbydup = window.slotbydup || []).push({
id: '4740881',
container: s,
size: '200,200',
display: 'inlay-fix'Matlab_学霸学习网
Matlab 软件目录1、Matlab 帮助的使用 ........................................................................ 3 2、Matlab 数据输入与类型 ................................................................ 4 3、Matlab 中的 M 文件及程序调试 ................................................ 12 4、Matlab 绘图命令 .......................................................................... 17 5、Matlab 在高等数学中的应用...................................................... 34 6、Matlab 在线性代数中的应用...................................................... 60 7、Matlab 数据处理 .......................................................................... 67 9、评价方法 ...................................................................................... 82 10、预测方法 .................................................................................... 97 11、蒙特卡洛方法 .......................................................................... 110 12、智能算法 .................................................................................. 122 13、分形 .......................................................................................... 129 14、Simulink 初步 .......................................................................... 134 15、Matlab 在概率统计中的应用.................................................. 147 参考文献 .......................................................................................... 1511 作为和 Mathematica、 Maple 并列的三大数学软件。 其强项就是其强大的矩阵计算以 及仿真能力。要知道 Matlab 的由来就是 Matrix + Laboratory = Matlab，所以这个软件在 国内也被称作《矩阵实验室》 。每次 MathWorks 发布 Matlab 的同时也会发布仿真工具 Simulink。在欧美很多大公司在将产品投入实际使用之前都会进行仿真试验，他们所主 要使用的仿真软件就是 Simulink。 Matlab 提供了自己的编译器： 全面兼容 C++以及 Fortran 两大语言。所以 Matlab 是工程师，科研工作者手上最好的语言，最好的工具和环境。 Matlab 已经成为广大科研人员的最值得信赖的助手和朋友！ 目前 MATLAB 产品族可以用来进行： － 数值分析 － 数值和符号计算 － 工程与科学绘图 － 控制系统的设计与方针 － 数字图像处理 － 数字信号处理 － 通讯系统设计与仿真 － 财务与金融工程... Simulink 是基于 MATLAB 的框图设计环境，可以用来对各种动态系统进行建模、 分析和仿真，它的建模范围广泛，可以针对任何能够用数学来描述的系统进行建模，例 如航空航天动力学系统、卫星控制制导系统、通讯系统、船舶及汽车等等，其中了包括 连续、离散，条件执行，事件驱动，单速率、多速率和混杂系统等等。 Simulink 提供 了利用鼠标拖放的方法建立系统框图模型的图形界面，而且 Simulink 还提供了丰富的 功能块以及不同的专业模块集合，利用 Simulink 几乎可以做到不书写一行代码完成整 个动态系统的建模工作。2 1、Matlab 帮助的使用1.1 helphelp?????????????%帮助总览 help elfun ???%关于基本函数的帮助信息 help exp ?????%指数函数 exp 的详细信息1.2 lookfor 指令当要查找具有某种功能但又不知道准确名字的指令时， help 的能力就不够了， lookfor 可以根据用户提供的完整或不完整的关键词，去搜索出一组与之相关的指令。 lookfor integral ??%查找有关积分的指令 lookfor fourier ??%查找能进行傅里叶变换的指令1.3 超文本格式的帮助文件在 Matlab 中，关于一个函数的帮助信息可以用 doc 命令以超文本的方式给出，如?doc ??? doc doc ??? doc eig ????%eig 求矩阵的特征值和特征向量1.4 pdf 帮助文件可从 MathWorks 网站上下载有关的 pdf 帮助文件。 网站地址：/3 2、Matlab 数据输入与类型2.1 Matlab 中的变量MATLAB 程序中的基本数据单元称为阵列(Array)，是一个分为行与列的数据集合。 变量被看作是只有一行一列的阵列。MATLAB 语言不需要对变量进行事先声明，也不 需要指定变量类型， 它会自动根据所赋予变量的值或对变量所进行的操作来确定变量的 类型。其命名规则为： （1）变量名的大小写是敏感的。 （2）变量的第一个字符必须为英文字母，而且不能超过 31 个字符。 （3）变量名可以包含下划线、数字，但不能为空格符、标点。 （4）命名变量时可以取一个容易记忆并且能表达出其含义的名称，如汇率，可以 定义为 exchange_rate。 对于变量作用域，默认情况是局部变量，使用 global 定义全局变量，而且全局变 量常用大写的英文字母表示。 MATLAB 预定义的变量如下表所示： ans eps pi inf NaN i 或 j nargin nargout realmax realmin flops 预设的计算结果的变量名 MATLAB 定义的正的极小值=2.2204e-16 内建的π 值 ∞值，无限大 无法定义一个数目 虚数单位 i=j=√-1 函数输入参数个数 函数输出参数个数 最大的正实数 2 最小的正实数 2 浮点运算次数41023-1022 注 1：在定义变量时要尽量与避免与这些名字相同，以免改变它们的值，如果已经 改变，可以通过 clear 复这些值。 注 2：数字的输入输出格式。所有数据均按 IEEE 浮点标准的长型格式存储。输入 格式沿用了 C 语言的风格和规则；输出格式使用 format 数据格式 命令控制，只影响 在屏幕上的显示结果，不影响内部的存储和运算。 变量名 来恢复它的初始值，也可以通过重新启动 MATLAB 恢2.2 向量及其运算1、向量的生成 ①命令窗口直接输入，使用[ ]，元素之间用空格、逗号或者分号隔开。 ②使用冒号表达式，基本形式为 x=x0:step:xn，其中 xn 为尾元素数值限，而不一定 是尾元素的值。当 step=1 时可省略步长。 ③生成线性等分向量，使用 linspace 函数。Y=linspace(x1,x2,n) ④生成对数等分向量，使用 logspace 函数。Y=logspace(x1,x2,n) 2、向量的基本运算5 数加（减） 、数乘、点积(dot 函数)、叉积(cross 函数)、混合积 dot(a,cross(b,c))2.3 矩阵及其运算1、简单矩阵的输入 （1）要直接输入矩阵时，矩阵一行中的元素用空格或逗号分隔；矩阵行与行之间 用分号“;”隔离，整个矩阵放在方括号“[ ]”里。 A=[1,2,3;4,5,6;7,8,9] ??? 说明：指令执行后，矩阵 A 被保存在 Matlab 的工作空间中，以备后用。如果用户 不用 clear 指令清除它，或对它进行重新赋值，那么该矩阵会一直保存在工作空间中， 直到本次指令窗关闭为止。 （2）矩阵的分行输入，此时回车键作为分行标志， A=[1,2,3 4,5,6 7,8,9] （3）使用 M 文件创建大矩阵，当矩阵维数非常大时，可以创建 m 文件，在 m 文 件中输入数据或者导出数据文件。 2、矩阵的基本运算 ①矩阵的四则运算。其中乘法运算要注意相乘的双方有相邻公共维，除法分为左除 “\” (A\B=inv(A)*B)和右除”/”（ A/B=A*inv (B)）(需要计算逆矩阵) ②矩阵的逆运算。inv 函数。 ③矩阵的幂运算。^。 ④矩阵的指数运算。exp(返回每个元素的指数值)，expm([V,D] = EIG(X) and EXPM(X) = V*diag(exp(diag(D)))/V)，expm1(exp(x)-1) ⑤矩阵的对数运算。logm ⑥矩阵的特征值函数。eig 和 eigs(适合于大型稀疏方阵) ⑦矩阵的奇异值函数。svd（[U,S,V] = SVD(X)，X = U*S*V'）和 svds ⑧矩阵的条件数函数。cond（矩阵 A 的条件数等于 A 的范数与 A 的逆的范数的 乘积，c = cond(A,p)等价于 norm(A,p) * norm(inv(A),p)） ，condest(1 范数的条件数的估计 值)， rcond ⑨特征值的条件数函数。codeig（[V,D,s] = condeig(A) 等价于[V,D] = eig(A); s6 =condeig(A);） ⑩范数函数。norm（1-范数：即列范数，矩阵的各列绝对值之和的最大值；2-范数： 所有元素的平方和开根号(默认)；无穷范数：即行范数，矩阵各行的绝对值之和的最大 值） ， normest (矩阵的 2 范数的估计值) 其他还有秩函数 rank， 迹函数 trace， 零空间函数 null （又称为核空间， X=null(A)， 则 A*X=0，X?*X=I） ，正交空间函数 orth(B = orth(A),则 B'*B = eye(rank(A)))，伪逆函数 pinv 等。 3、特殊向量和特殊矩阵 （1）特殊向量 t=[0:0.1:10] %产生从 0 到 10 的行向量，元素之间间隔为 0.1 t=linspace(n1,n2,n) %产生 n1 和 n2 之间线性均匀分布的 n 个数 (缺省 n 时,产生 100 个数) t=logspace(n1,n2,n) (缺省 n 时,产生 50 个数) %在和之间按照对数距离等间距产生 n 个数。 （2）特殊矩阵 i)单位矩阵 eye(m), eye(m,n) 可得到一个可允许的最大单位矩阵而其余处补 0, eye(size(a)) 可以得到与矩阵 a 同样大小的单位矩阵。 ii) 所有元素为 1 的矩阵 ones(n)，ones(size(a))，ones(m， n)。 iii）所有元素为 0 的矩阵 zeros(n)，zeros(m,n)。 iv）空矩阵是 一个特殊矩阵，这在线性代数中是不存在的。例如 q=[ ] 矩阵 q 在工作空间之中，但它的大小为零。通过空矩阵的办法可以删除矩阵的行与 列。例如 a(:,3)=[] 表示删除矩阵 a 的第 3 列。 v）随机数矩阵7 rand(m,n) 产生 m× n 矩阵，其中的元素是服从[0,1]上均匀分布的随机数。 randint(m,n,[min,max]) 产生 m× n 矩阵，其中的元素是[min,max]上的随机整数。 normrnd(mu,sigma,m,n)产生 m× n 矩阵，其中的元素是服从均值为 mu，标准差 为 sigma 的正态分布的随机数。 exprnd(mu,m,n) 产生 m× n 矩阵，其中的元素是服从均值为 mu 的指数分布的随 机 数。 poissrnd(mu,m,n) 产生 m× n 矩阵，其中的元素是服从均值为 mu 的泊松 （Poisson） 分布的随机数。 unifrnd(a,b,m,n) 产生 m× n 矩阵， 其中的元素是服从区间[a,b]上均匀分布的随机数。 r = mvnrnd(MU,SIGMA,cases) 产生 cases 对均值向量为 MU，协方差阵为 SIGMA 的多维正态分布的随机数。 vi）随机置换 randperm(n)产生 1 到 n 的一个随机全排列。 perms([1:n])产生 1 到 n 的所有全排列。 vii) 稀疏矩阵 稀疏矩阵是指矩阵中零元素很多，非零元素很少的矩阵。对于稀疏矩阵，只要存放 非零元素的行标、列标、非零元素的值即可，可以按如下方式存储 （非零元素的行地址，非零元素的列地址），非零元素的值。 在 Matlab 中无向图 和有向图邻接矩阵的使用上有很大差异。 对于有向图，只要写出邻接矩阵，直接 使用 Matlab 的命令 sparse 命令，就可以把 邻接矩阵转化为稀疏矩阵的表示方式。 对于无向图，由于邻接矩阵是对称阵，Matlab 中只需使用邻接矩阵的下三角元素， 即 Matlab 只存储邻接矩阵下三角元素中的非零元素。 稀疏矩阵只是一种存储格式。Matlab 中，普通矩阵使用 sparse 命令变成稀疏矩阵， 稀疏矩阵使用 full 命令变成普通矩阵。 例 1 a=zeros(5); a(1,[2,4])=[3,4]; a(3,[2:4])=[1 3 8]; a(5,[1,5])=[6,7]8 b=sparse(a) %普通矩阵转化成稀疏矩阵 c=full(b) %稀疏矩阵转化成普通矩阵 其他一些特殊矩阵如下表所示： 函数 compan gallery hadamard hankel hilb invhilb 功能 伴随阵 Higham 测试阵 Hardamard 矩阵 Hankel 矩阵 Hilbert 矩阵 反 Hilbert 矩阵 函数 magic rosser toeplitz pascal vander wilkinson 功能 魔方阵 经典对称特征值测试阵 Toeplitz 矩阵 Pascal 矩阵 范德蒙矩阵 Wilkinson’s 特征值测试矩阵4、矩阵的特殊操作 ①变维。 有两种方法，使用冒号（:）和使用函数 reshape 使用“:”表达式对两个矩阵进行变维操作，需要预先定义两个矩阵的维数（例如 A(:,3)=[]） ；reshape 有两种形式，分别为 reshape(X,M,N)和 reshape(X,M,N,P…) ②变向 主要函数如下表所示： 函数 fiplr fipud fipdim Rot90 ③矩阵的抽取 对角线元素抽取函数 diag(X,k)/diag(v,k)，抽取矩阵 X 的第 k 条对角线的元素向量/ 使得向量 v 为所得矩阵的第 k 条对角线元素。 上三角元 素抽取 tril(X,k)和下三角元素抽取 triu(X,k) ④扩展 两种方法：利用对矩阵标示块的赋值命令 X(m1:m2,n1:n2)=a 生成大矩阵， 其中 m2-m1+1 必须等于 a 的行维数，n2-n1+1 必须等于 a 的列维数，生成 m2 × n2 维的 功能 矩阵左右翻转 矩阵上下翻转 矩阵特定维翻转 矩阵反时针 90 翻转 函数 diag tril triu 功能 产生或提取对角阵 产生下三角 产生上三角9 矩阵 X；利用小矩阵组合生成大矩阵，要严格注意矩阵大小的匹配。 （示例： X(1:10,1:10)=zeros(10,10)，LX=[X,X;X,X]） Matlab 中冒号（： ）的使用方法小结： （1）用于生成向量，a:b，一般要求 a&b，否则生成空矩阵。 （2）如果 b-a 不是整数时，则向量的最后一位数为 n+a，其中 n=fix(b-a)(向零取整) （3）a:c:b 表示[a,a+c,…,a+n*c]，其中 n=fix((b-a)/c)，当 c&0 且 a&b 或者 c&0 且 a&b 时出现空矩阵。 （4）A(:)以一列的方式显示 A 中所有的元素 （5）b=A(i,:)表示将 A 中第 i 行存入 b 中 （6）b=A(:,j)表示将 A 中第 j 列存入 b 中 （7）b=A(j:k)表示将 A 中第 j 到第 k 个元素存入 b 中（Matlab 中矩阵按列存储） （8）b=A(:,c:d)表示将 A 中第 c 列到第 d 列存入 b 中，要求 c，d 不能超过 A 的列 数。 （9）当矩阵很大时，不知道矩阵的维数，可以使用 end 作为矩阵的最后一行或者 一列或者最后一个元素。例如 b=A(1:2,2:end)获取矩阵 A 右上角的元素。2.4 数组及其运算数组与矩阵在形式上完全一致，只是运算与矩阵不同。同型矩阵之间的运算通常称 为数组运算。 （矩阵的数组运算） 1、基本数组运算 ①四则运算。数组的乘除法是指两个同维数组间对应元素之间的乘除法，运算符 为”.*”,”./”和”.\”。数组与常数之间的运算可以加”.”，也可以不加。 ②幂运算。.^,对每个数组元素的幂运算。 ③指数运算 exp,对数运算 log 和开方运算 sqrt。 2、数组函数运算 只要把运算的数组带入到函数中就可以了，通用 形式为 funname(A) 3、数组的逻辑运算和关系运算 指令 & 含义 小于 函数名 lt10 &= & &= == ~= & | ~小于等于 大于 大于等于 等于 不等于 逻辑 与 逻辑 或 逻辑 非le gt ge eq ne and or not指令 xor any all isnan含义 不相同就取 1，否则取 0 只要有非 0 就取 1，否则取 0 全为 1 取 1，否则为 0 为数 NaN 取 1，否则为 0指令 isequal ismember isempty isletter含义 相等取 1，否则取 0 两个矩阵是属于关系取 1， 否则取 0 矩阵为空取 1，否则取 0 是字母取 1， 否则取 0 （可以是字符 串）isinf isfinite ischar find isnumeric为数 inf 取 1，否则为 0 有限大小元素取 1，否则为 0 是字符串取 1，否则为 0 寻找非零元素坐标 判断数值矩阵isstudent isprime isreal isspace islogical学生版取 1 质数取 1，否则取 0 实数取 1，否则取 0 空格位置取 1，否则取 0 判断逻辑数组11 3、Matlab 中的 M 文件及程序调试M 文件分为两种，一种是脚本文件，由一系列 Matlab 的命令组成，可以直接运行； 一种是函数文件，必须由其他 M 文件或者在命令行窗口中调用执行。函数文件具有一 定的通用性，并且可以进行递归调用。3.1、Matlab 中流程控制语句（1）if 语句 有三种形式：if(表达式) 语句组 A; end if(表达式) 语句组 A，else 语句组 B，end if(表达式 1) 语句组 A，elseif(表达式 2) 语句组 B, else 语句组 C，end （2）循环语句 while（表达式）语句组 A , end for k=初值：增量：终值 语句组 A，end 循环语句可以结合 break 命令来控制程序的执行顺序。 （3）switch 语句 switch 表达式(标量或者字符串) case 值 1 语句组 A case 值 2 语句组 B … otherwise 语句组 N end 用法示例： 示例一：输入数 n，判断其奇偶性 解 Matlab 程序如下： n=input(?n=?); if rem(n,2)==0 %求余数，与 mod 函数的区别，当同符号时两者相同，否则不同 A=?even?; else A=?odd?,12 end 示例二：求 Matlab 中的最大实数 解：Matlab 程序如下： x=1; while x~=inf, x1=x; x=2*x; %为了获取更接近 realmax 的值，可以改为 x=1.01*x end x1 示例三：求 Matlab 的相对精度 解：matlab 程序如下： y=1; while 1+y&1 y1=y; y=y/2; end y13.2 脚本文件主要特征如下： (1)一般由 clc,clear 命令开始，清除掉屏幕和工作空间中原有的变量和图形，以避免 其他已执行的程序残留数据对本程序的影响。程序中应该添加有注释。 (2)接下来是程序的主体，如果文件中有全局变量，则需要使用 Global 在程序的起 始部分注明。为了提高程序的可读性，注意语句的缩进。 (3)整个程序应按照 Matlab 标识符的要求起文件名，扩展名为 m。 示例：列出求素数的程序 解 Matlab 程序如下： clear,clc N=input('N='),x=2:N; %列出从 2 到 N 的全部自然数 for u=2:sqrt(N) %依次列出除数 n=find(rem(x,u)==0&x~=u);%找到能被 u 整除而不等于 u 的数序号 x(n)=[]; %剔除该数 end x13 3.3 函数文件函数文件与脚本文件的区别： (1)由 function 开头，后跟的函数名必须与文件名相同； (2)有输入变量和输出变量，可进行变量传递； (3)除非使用 Global 声明，程序中的变量均为局部变量，运行后不保存在工作空间 中。 示例 1：使用函数文件写出求素数的程序 解：Matlab 程序如下： function y=qiuprime(N) %求出 1 到 N 之间的素数，返回素数组成的向量 x=2:N; %列出从 2 到 N 的全部自然数 for u=2:sqrt(N) %依次列出除数 n=find(rem(x,u)==0&x~=u);%找到能被 u 整除而不等于 u 的数序号 x(n)=[]; %剔除该数 end y=x; 示例 2：写一个递归函数求 n! 解： function f=factor(n) if n&=1 f=1; else f=factor(n-1)*n; end end3.4 函数参数的可调性Matlab 中有两个永久变量 nargin 和 nargout 分别记录调用该函数时的输入实参和输 出实参的个数。只要在函数中包含这两个变量，就可以准确的知道该函数文件被调用时 的输入输出参数个数，从而决定函数如何进行处理。其他类似的变量还有 varargout， varargin。14 示例 1： function fout=examp(a,b,c) if nargin==1 fout=a; elseif nargin==2 fout=a+b; elseif nargin==3 fout=a+b+c; end 示例 2： function [x,y,z]=examp(a,b,c) if nargout==1 x=a; elseif nargout==2 x=a; y=b; elseif nargout==3 x=a; y=b; z=c; end end 示例 3： function vartest(argA, argB, varargin) optargin = size(varargin,2); %可选输入 stdargin = nargin - %标准输入 fprintf('Number of inputs = %d\n', nargin) fprintf(' Inputs from individual arguments(%d):\n', ... stdargin) if stdargin &= 1 fprintf(' %d\n', argA) end if stdargin == 2 fprintf(' %d\n', argB) end fprintf(' Inputs packaged in varargin(%d):\n', optargin) for k= 1 : size(varargin,2) fprintf(' %d\n', varargin{k}) end end15 3.5 程序调试M 文件编辑器中提供了强大的程序调试功能，使用方式如同 VC 程序中的调试器。 在 M 文件编辑器中的工具栏中有如下程序调试图标：这些图标的使用和菜单 Debug 中的一些命令相同。16 4、Matlab 绘图命令4.1 二维绘图命二维绘图的基本命令有 plot，loglog，semilogx，semilogy 和 polar。它们的使用方法 基本相同，其不同特点是在不同的坐标中绘制图形。plot 命令使用线性坐标空间绘制图 形；loglog 命令在两个对数坐标空间中绘制图形；而 semilogx(或 semilogy)命令使用 x 轴（或 y 轴）为对数刻度，另外一个轴为线性刻度的坐标空间绘制图形；polar 使用极 坐标空间绘制图形。 二维绘图命令 plot 为了适应各种绘图需要，提供了用于控制线色、数据点和线型的 3 组基本参数。它的使用格式如下： plot(x,y,‘color_point_linestyle’) 该命令是绘制 y 对应 x 的轨迹的命令。y 与 x 均为向量，且具有相同的元素个数。 用字符串，color_point_linestyle? 完成对上面 3 个参数的设置。 线色：r-red，g-green，b-blue，w-white，k-black，i-invisible，y-yellow，c-cyan （青色） ，m-紫色。 数据点：.（圆点） ，o（小圆圈） ，x（叉号） ，+（加号） ，*（星号） ，s （square 方形） ， h（hexagram 六角星） ，d（diamond 菱型） ，p（pentagram 五角星） ，v（下三角） ，^（上 三角） ，&（左三角） ，&（右三角） 。 线型：-（实线） ，--（虚线） ，-.（点画线） ，:（点线） 。 当 plot(x,y)中的 x 和 y 均为 m× n 矩阵时，plot 命令将绘制 n 条曲线。 plot(t,[x1,x2,x3])在同一坐标轴内同时绘制三条曲线。 如果多重曲线 对应不同的 x 轴向量绘制，可使用命令 plot(t1,x1,t2,x2,t3,x3) 式中 x1 对应 t1，x2 对应 t2 等等。在这种情况下，t1，t2 和 t3 可以具有不同的元素 个数，但要求 x1，x2 和 x3 必须分别与 t1，t2 和 t3 具有相同的元素个数。 subplot 命令使得在一个屏幕上可以分开显示 n 个不同坐标系，且可分别在每一个 坐标系中绘制曲线。其命令格式如下 subplot(r,c,p)17 该命令将屏幕分成 r× c 个子窗口，而 p 表示激活第 p 个子窗口。窗口的排号是从 左 到右，自上而下。 下面对几个特殊的坐标系统 进行简要介绍： ①对数坐标曲线，主要有 semilogx,semilogy 和 loglog，前两个分别以 x 坐标和 y 坐标为对数坐标，后一个是双对数坐标。 例如： x=1:0.1*pi:2* y=sin(x); semilogx(x,y,'-*') x = 0:.1:10; semilogy(x,10.^x) x = logspace(-1,2); loglog(x,exp(x),'-s') grid on ②双纵坐标（双 y 轴坐标系）函数 plotyy，调用形式为： plotyy(X1,Y1,X2,Y2) plotyy(X1,Y1,X2,Y2,fun) fun 可以是 plot、semilogx、semilogy 或 loglog plotyy(X1,Y1,X2,Y2,fun1,fun2) fun1 绘制(X1,Y1)，fun2 绘制(X2,Y2) ③极坐标绘图函数 polar 调用形式为：polar(theta,rho) 例 1：绘制极坐标下的平面曲线 ? ? a cos(b ? n? ) ，并讨论参数 a,b,n 对曲线的影响。 解 Matlab 程序如下： N=100; theta=linspace(0,2*pi,N); for i=1:2 a(i)=input('a='); b(i)=input('b='); n(i)=input('n='); rho(i,:)=a(i)*cos(b(i)+n(i)*theta); subplot(1,2,i),polar(theta,rho(i,:)); end 几个比较漂亮的极坐标系下的图 形： (1)蝴蝶图案18 t=0:0.01:36; f=exp(cos(t-pi/2))-2*cos(4*(t-pi/2))+sin((t-pi/2)/12).^5; polar(t,f,'r') (2)枫叶图案 t=-pi/2:0.05:1.5* f=100./(100+(t-pi/2).^8)*2-sin(7*t)-cos(30*t)/2; polar(t,f,'r') 在图形绘制完毕后， 执行如 下命令可以再在图中加入标题、 标号、 说明和分格线等。 这些命令有 title，xlabel，ylabel，text，gtext，legend 等。它们的命令格式如下 title(“My Title”)，xlabel(“My X-axis Label ”)，ylabel(“My Y-axis Label”)， xlabel({'first line';'second line'})， text(x, y,'Text for annotation')， gtext('Text for annotation')， grid on gtext 命令是使用鼠标器定位的文字注释命令。当输入命令后，可以在屏幕上得到 一个光标， 然后使用鼠标器控制它的位置。 按鼠标器的左键， 即可确定文字设定的位置。 hold on 是图形保持命令，可以把当前图形保持在屏幕上不变，同时在这个坐标系 内绘制另外一个图形。 hold on 命令是一个交替转换命令， 即执行一次， 转变一个状态 （相 当于 hold on、hold off） 。 还可以设置坐标轴的范围，使用命令 axis，其格式为： axis([xmin xmax ymin ymax]) axis([xmin xmax ymin ymax zmin zmax cmin cmax]) 示例： x = 0:.01:pi/2; plot(x,tan(x),'-ro') 图形效果如下所示：默认情况下，横纵坐标轴的范围是根据函数自变量和因变量的值自动变化的，有时 效果不好，此时需要设定横纵坐标轴的范围：axis([0 pi/2 0 5])，效果如上右图所示。19 解 画图的 Matlab 程序如下，画出的图见图 1。 clc, clear x=-2*pi:0.1:2* y1=sin(x); y2=sin(x+pi/3)+2; y3=cos(x); plot(x,y1,'.-'); hold on %图形保持命令 plot(x,y2,'*-'); plot(x,y3,'-o'); h=legend('sin($x$)','sin($x+\frac{\pi}{3}$)','cos($x$)') %latex 格式显示 set(h,'Interpreter','latex') %设置 Interpreter 的属性值为 latex，可以使用数学公式 xlabel('$x$','Interpreter','latex') %latex 格式显示 ylabel('$y$','Interpreter','latex') %latex 格式显示图1解 画图的 Matlab 程序如下，画出的图形结果见图 2。 clc, clear x=-2*pi:0.1:2* y1=sin(x); y2=sin(x+pi/3)+2; y3=cos(x);20 subplot(3,1,1), plot(x,y1,'.-'), title('sin($x$)','Interpreter','latex') subplot(312), plot(x,y2,'*-') title('sin($x+\frac{\pi}{3}$)','Interpreter','latex') ylabel('$y$','Interpreter','latex') %latex 格式显示 subplot(313), plot(x,y3,'-o') title('cos($x$)','Interpreter','latex') %latex 格式显示 xlabel('$x$','Interpreter','latex') %latex 格式显示4.2 复数绘图plot(z) ：z 为复数时相当于 plot(real(z),imag(z))；如果是双变量如 plot(t,z)，则 z 中 的虚数部分将丢弃。在复平面中绘出多条曲线，必须使用 hold on 命令，或者把多条曲 线的实部和虚部明确的写出。即 plot(real(z1),imag(z1),real(z2),imag(z2)) Matlab 中专门用来绘制复变量函数图形的相关命令是 cplxmap， cplxgrid， cplxroot， 格式如下： z = cplxgrid(m) %产生(m+1)*(2m+1)的极坐标下的复数数据网格 cplxmap(z,f(z),(optional bound)) %绘制复变函数的图形，以 xy 平面为自变量所在的 复平面，以 z 轴表示复变函数的实部，颜色表示复变函数的虚部。 cplxroot：画复数的 n 次函数曲面21 cplxroot(n) %画复数 n 次根的函数曲面,复数为最大半径为 1 的圆面 cplxroot(n,m) %画复数 n 次根的函数曲面,复数为最大半径为 1 的圆面,为（m+1）* （2m+1）的方阵 例 1：绘制 z=exp((-0.1+i)*t)的复数图形 解 Matlab 程序如下： t=0:0.1:15; z=exp((-0.1+i)*t); subplot(2,2,1) plot(z), % pause 可以暂停程序的执行， 通过按任意键可使程序继续进行， 也 可以在 pause(n)中设置时间使执行结果出现动态效果。 title('复数绘图 plot(z)'); subplot(2,2,2) plot(t,z),pause title('plot(t,z)') subplot(2,2,3) polar(angle(z),abs(z)); title('polar(angle(z),abs(z))') subplot(2,2,4) semilogx(t,z); title('semilogx(t,z)') 注：Matlab 中使用函数 C=complex(A,B)构造复数。 例 2：绘图 y=cos(x+i)的图形 解法一：x=-pi:0.1: fun=@(x)cos(x+i); %如果直接使用 plot(x,y)，会忽略虚部。 plot(fun(x)) 解法二：plot(cos(x+i)); 例 3：绘制 f=z^4 的图形，其中 z 为复数 解：z=cplxgrid(30); cplxmap(z,z.^4)4.3 显函数，符号函数或隐函数的绘图fplot(fun,lims)绘制由字符串 fun 指定函数名的函数在 x 轴区间为 lims=[xmin, xmax] 的函数图。若 lims=[xmin,xmax,ymin,ymax]，则 y 轴也被限制。22 解 （1）首先用 M 文件 fun1.m 定义函数 f(x)如下 function y=Afun1(x); if x&1 y=x+1; else y=1+1./x; end 在 matlab 命令窗口输入 fplot('Afun1',[-3,3]) 就可画出 函数 f (x) 的图形。 （2）可以使用匿名函数，编写程序如下 fun2=@(x) (x+1)*(x&1)+(1+1/x)*(x&=1); fplot(fun2,[-3,3]) ezplot(f)绘制符号函数或者隐函数 f(x)的图形，x 轴的近似范围为[2?,?2??。 ezplot(f,[xmin,xmax])使用输入参数来代替默认横坐标范围[2?,?2?? 。 ezplot 函数的其 他格式有： ? ezplot(fun2)绘制 fun2(x,y)=0 的隐函数曲线，默认 x,y 的范围是[2 ?,?2??? ezplot(fun2,[xmin,xmax,ymin,ymax]) 示例： function z = myfun(x,y,k) %建立 M 文件 z = x.^k - y.^k - 1; %在命令窗口中输入：ezplot(@(x,y)myfun(x,y,2));解 ezplot('cot(x)')解 ezplot('x^2+y^2/4=1') 例 6-1：同一坐标系下绘制曲线 x^2+y^2=1 和 x^2-y^2=1 的曲线 解：ezplot('x^2-y^2=1');23 ezplot('x^2+y^2=1'); colormap([0,0,1]);4.4 三维图形在实际工程计算中，最常用的三维绘图是三维曲线图、三维网格图和三维曲面图 3 种基本类型。 与此对应， Matlab 也提供了一些三维基本绘图命令， 如三维曲线命令 plot3， 三维网格图命令 mesh 和三维表面图命令 surf。 1. 三维曲线 plot3(x,y,z)通过描点连线画出曲线， 这里 x,y,z 都是 n 维向量， 分别表示该曲线上点 集的横坐标、纵坐标、竖坐标。 例7 t=0:0.05:20* x=sin(t); y=cos(t); z=t.*sin(t).*cos(t); plot3(x,y,z), title('Line in 3-D Space') xlabel('X'), ylabel('Y'), zlabel('Z'), grid on 2. 三维网格图 命令 mesh(x,y,z)画网格曲面。这里 x,y,z 是三个同维数的数据矩阵，分别表示数据 点的横坐标、纵坐标、竖坐标，命令 mesh(x,y,z)将该数据点在空间中描出，并连成网格。解 x=-3:0.1:3;y=-5:0.1:5; x1=ones(size(y'))*x;y1=y'*ones(size(x)); [x2,y2]=meshgrid(x,y); %meshgrid 函数生成 2D/3D 网格矩阵，用于分割空间，在绘 制 3D 网格图或者表面图时都需要调用该函数对生成绘图时所用的数据。 z1=(sin(x1.*y1)+eps)./(x1.*y1+eps); z2=(sin(x2.*y2)+eps)./(x2.*y2+eps);24 subplot(1,2,1),mesh(x1,y1,z1) subplot(1,2,2),mesh(x2,y2,z2) 例 8-1：绘制两个空间相交的平面 z=x+2y-1 和 z=x-y+2 解 ： x=-3:0.1:3;y=-3:0.1:3; [X,Y]=meshgrid(x,y); mesh(X,Y,X+2*Y-1) hold on mesh(X,Y,X-Y+2) 3. 表面图 命令 surf(x,y,z)画三维表面图，这里 x,y,z 是三个同维数的数据矩阵，分别表示数据 点的横坐标、纵坐标、竖坐标。解 [x,y]=meshgrid([-3:0.2:3]); z=(sin(x.*y)+eps)./(x.*y+eps); surf(x,y,z) 4. 旋转曲面 方法一：使用命令函数 cylinder [X,Y,Z] = cylinder(r) 这里的 r 表示构成旋转曲面的曲线。解 Matlab 程序如下。 x=0:10:600; [X,Y,Z]=cylinder(30*exp(-x/400).*sin((x+25*pi)/100)+130); surf(X,Y,Z) 方法二：将旋转曲面用参数方程表示。25 解 因为这里的函数是隐函数，化成显函数后有两支，必须使用参数方程，旋转面 的参数方程为画图的 Matlab 程序如下： alpha=[0:0.1:2*pi]'; beta=0:0.1:2* x=4*cos(alpha)*ones(size(beta)); y=(5+4*sin(alpha))*cos(beta); z=(5+4*sin(alpha))*sin(beta); surf(x,y,z) 或者利用绘制三维隐函数曲面图形的命令 ezsurf 或者 ezmesh，其命令格式如下： ezmesh(fun) ezmesh(fun,domain) domain 要求是一个向量，指定 x,y 轴的范围 ezmesh(funx,funy,funz) %参数方程绘图 funx(s,t), funy(s,t), and funz(s,t) over the square: -2π & s & 2π, -2π & t & 2π. ezsurf(fun) ezsurf(fun,domain) ezsurf(funx,funy,funz) 画图的 Matlab 程序也可以写成 x=@(alpha,beta) 4*cos(alpha); y=@(alpha,beta) (5+4*sin(alpha))*cos(beta); z=@(alpha,beta) (5+4*sin(alpha))*sin(beta); ezsurf(x,y,z) 对于形如 f(x,y,z)=0 的三维隐函数的图像没有现成的函数可以画，但可以利用 isosurface 函数绘制三角网格图。例如 f=@(x,y,z)x.^2+y.^2+z.^2-10;%定义函数 f=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,0);%用 isosurface 得到函数 f=0 图形的点和面26 patch('faces',p,'vertices',v,'facevertexcdata',jet(size(v,1)),'facecolor','w','edgecolor','flat'); %用 patch 绘制三角网格图并设定色彩 view(3); axis equal 5、绘制柱面 柱面平行于某个坐标轴，方程中不出现某个坐标轴的变量，方程表示为 F(x,y)=0 或 者 F(x,z)=0 或者 F(y,z)=0 示例一：画出方程 z ? 2 ? x2 表示的柱面 解：方程中的 x 是自变量矩阵，z 是因变量，则另一个自变量矩阵为 y，自变量平 面是 xoy 面，x 轴是真正的自变量，y 轴是柱面方向。Matlab 程序如下： u=linspace(-5,5,10)';%设定参数列向量 u v=linspace(-5,10,10); %设定参数行向量 v X=u*ones(size(v)); %构成自变量矩阵 X Y=ones(size(u))*v; %构成自变量矩阵 Y Z=2-X.^2; %求因变量 Z mesh(X,Y,Z) 示例二：画出方程 ? x 2 ? y2 ? 1 表示的柱面 4解：方程整理为显示函数： y ? ?2 1 ? x2 ，分别绘制正负两个分量。 clear,clc u=linspace(-5,5,10)';%设定参数列向量 u v=linspace(-5,10,10); %设定参数行向量 v X1=u*ones(size(v)); %构成自变量矩阵 X1 Y1=2*sqrt(1+X1.^2); %求因变量 Y1 Y2=-2*sqrt(1+X1.^2); Z1=ones(size(u))*v; %求因变量 Z1 mesh(X1,Y1,Z1),hold on mesh(X1,Y2,Z1) 6. 其它二次曲面 对于旋转面，如果母线的方程可以表示成关于旋转轴变量的显式函数，则可以直接 使用 Matlab 工具箱中的命令 cylinder，否则必须把旋转面化成参数方程，然后使用 ezmesh 或 ezsurf 命令绘图。对于其它的二次曲面，如果可以写成显函数直接使用命令 ezmesh 或 ezsurf，否则必须先化成参数方程。还有一些特殊的绘制函数如27 cylinder,ellipsoid 等命令。解： （1）x=@(s,t) 3*sec(s)*cos(t); %化为参数方程 y=@(s,t) 3*sec(s)*sin(t); z=@(s,t) 2*tan(s); ezmesh(x,y,z) （2）x=@(s,t) 3*sec(s); y=@(s,t) 2*tan(s)*cos(t); z=@(s,t) 2*tan(s)*sin(t); ezmesh(x,y,z) （3）ezsurf(@(y,z) y.^2,50) %直接调用 ezsurf （4）x=@(s,t) 3*tan(s)*cos(t); y=@(s,t) 2*tan(s)*sin(t); z=@(s,t) tan(s); ezsurf(x,y,z) （5）ellipsoid(0,0,0,3,2,sqrt(6)) %专门绘制椭球面 （6）ezsurf(@(x,y)x*y) （7）x=@(s,t) 3*cos(s);28 y=@(s,t) 2*sin(s); z=@(s,t) ezmesh(x,y,z) 7、绘制空间两曲面的交线 示例：绘制由水平截面与方程 z ? x2 ? 2 y 2 构成的马鞍面形成的交线，并讨论等高线 和方向导数（梯度）的意义 解： 水平平面与曲面的交线就是等高线， 在 Matlab 中绘制等高线有两个命令， contour 和 contour3，前者把等高线画在 xoy 平面上，后者把等高线画在一定高度的平面上，使 之成为立体的，与所在曲面对应。 Matlab 程序如下： clc,clear [x,y]=meshgrid(-10:.2:10);%确定计算和绘图的定义域网格 z1=(x.^2-2*y.^2)+ %第一个曲面方程 a=input('a=(-50&a&50)'); z2=a*ones(size(x)); %水平面方程 z2=a，z2 必须与 x，y 具有相同的维数， subplot(1,3,1),mesh(x,y,z1);mesh(x,y,z2);%分别画出两个曲面 v=[-10 10 -10 10 -100 100];axis(v), %确定第一个分图的坐标系 colormap(gray); r0=abs(z1-z2)&=1; %求两曲面 z 坐标只差小于 0.5 的网格 zz=r0.*z2;yy=r0.*y;xx=r0.*x; %求这些网格上的坐标值，即交线坐标值 subplot(1,3,2),plot3(xx,yy,zz,'x'); %画出这些点 axis(v), %使第二个分图取第一个分图的坐标系 pause,subplot(1,3,3); contour3(x,y,z1,20); %用等高线命令求出 20 条不同高度的交线。 等高线与方向导数和梯度的概念密切相关，函数 z1 在每一点的梯度与该处的等高 线垂直，也就是指向最陡的方向。Matlab 中求梯度用 gradient 函数，quiver 函数画出梯 度向量，这两个函数要求在给定的点阵上求梯度和画梯度向量。Matlab 程序如下： [x,y]=meshgrid(-10:2:10); z1=(x.^2-2*y.^2)+ contour(x,y,z1,20); [px,py]=gradient(z1,2,2); %以步长为 2 求 z1 的梯度的 x，y 分量 quiver(x,y,px,py); %绘制梯度向量 8、离散点的图形绘制相关函数 griddata，TripScatterInterp，scatter（scatter3） griddata 用来对离散数据进行曲面拟合 (1)ZI = griddata(x,y,z,XI,YI) 用二元函数 z=f(x,y)的曲面拟合有不规则的数据向量29 x,y,z。griddata 将返回曲面 z 在点（XI,YI）处的插值。曲面总是经过这些数据点（x,y,z） 的。输入参量（XI,YI）通常是规则的格点（像用命令 meshgrid 生成的一样） 。XI 可以 是一行向量，这时 XI 指定一有常数列向量的矩阵。类似地，YI 可以是一列向量，它 指定一有常数行向量的矩阵。 (2)[XI,YI,ZI] = griddata(x,y,z,xi,yi) 返回的矩阵 ZI 含义同上， 同时， 返回的矩阵 XI,YI 是由行向量 xi 与列向量 yi 用命令 meshgrid 生成的。 (3)[XI,YI,ZI] = griddata(.......,method) 用指定的算法 method 计算： ?linear ?：基于三角形的线性插值（缺省算法） ； ?cubic?： 基于三角形的三次插值； ?nearest?：最邻近插值法； 示例： x = rand(100,1)*4-2; y = rand(100,1)*4-2; z = x.*exp(-x.^2-y.^2); ti = -2:.25:2; [xi,yi] = meshgrid(ti,ti); zi = griddata(x,y,z,xi,yi); mesh(xi,yi,zi), hold on, plot3(x,y,z,'o'), hold off 在新版本的 Matlab 中，该函数逐渐被 TriScatteredInterp 函数所代替， TriScatteredInterp 的示例如下： x = rand(100,1)*4-2; y = rand(100,1)*4-2; z = x.*exp(-x.^2-y.^2); F = TriScatteredInterp(x,y,z); ti = -2:.25:2; [qx,qy] = meshgrid(ti,ti); qz = F(qx,qy);%计算指定位置的插值 mesh(qx,qy,qz); plot3(x,y,z,'o'); Matlab 中绘制散点图的命令是 scatter（scatter3） ，命令格式如下： scatter(X,Y,S,C) %X，Y 对应散点值，S 代表标记大小，C 代表颜色值，它都可以是 向量。 scatter3(X,Y,Z,S,C) 示例：30 A=[1.486,3.059,0.1;2.121,4.041,0.1;2.570,3.959,0.1;3.439,4.396,0.1; 4.505,3.012,0.1;3.402,1.604,0.1;2.570,2.065,0.1;2.150,1.970,0.1; 1.794,3.059,0.2;2.121,3.615,0.2;2.570,3.473,0.2;3.421,4.160,0.2; 4.271,3.036,0.2;3.411,1.876,0.2;2.561,2.562,0.2;2.179,2.420,0.2; 2.757,3.024,0.3;3.439,3.970,0.3;4.084,3.036,0.3;3.402,2.077,0.3; 2.879,3.036,0.4;3.421,3.793,0.4;3.953,3.036,0.4;3.402,2.219,0.4; 3.000,3.047,0.5;3.430,3.639,0.5;3.822,3.012,0.5;3.411,2.385,0.5; 3.103,3.012,0.6;3.430,3.462,0.6;3.710,3.036,0.6;3.402,2.562,0.6; 3.224,3.047,0.7;3.411,3.260,0.7;3.542,3.024,0.7;3.393,2.763,0.7]; x=A(:,1);y=A(:,2);z=A(:,3); scatter(x,y,5,z)%散点图4.5 动态可视化图形Matlab 中的动画命令，moviein，getframe 和 movie，用 getframe 把 Matlab 产生的 图形存储下来，每个图形成一个很大的列向量；再用 N 行这样的列保存 N 幅图，成为 一个大矩阵；movie 命令把他们连接起来重放，产生动画效果；moviein 用来预留存储空 间以加快运行的速度。有些情况也需要借助 pause 命令和循环语句来实现点轨迹的动态 演示。绘制彗星图的命令 comet 或者 comet3 也可以实现点的运动轨迹，不过这种方式 速度较快。 示例 1： axis equal, %把坐标设成相等比例 M=moviein(16);%为变量 M 预留 16 幅图的存储空间 for j=1:16 plot(fft(eye(j+16))); M(:,j)= end movie(M,30) %以每秒 30 帧的速度播放 M 中的图形 示例 2：使用 comet/comet3 绘制彗星图 （1） t = 0:.01:2* %返回当前坐标轴下的一帧图像31 x = cos(2*t).*(cos(t).^2); y = sin(2*t).*(sin(t).^2); comet(x,y); （2） n=10; t=n*pi*(0:0.0005:1); x=sin(t); y=cos(t); plot(x,y,'g'); comet(x,y,0.01); 示例 3：极坐标下一个点沿着图形移动的轨迹 a=0:.01:2* b=3; polar(a,b*(1-cos(a)),'*r'); for a=0:0.05:2*pi cla(findobj(gca,'color','r')); %findobj 函数用来使用特定的属性获取图形句柄， cla 用来清理由函数句柄指定的图形所在的坐标轴。 h2=polar(a,b*(1-cos(a)),'b'); set(h2,'linestyle','.','markersize',30); %设置句柄指定函数图形的属性 pause(0.001) end 示例 4： 普通二维坐标下一个点沿着曲线移动的轨迹， 程序与上面的结果基本相同。 x=-2*pi:0.1:2* y=sin(x); plot(x,y,'b-'); hold on for m=-2*pi:0.05:2*pi cla(findobj(gca,'color','b')); %cla 清除当前坐标轴上面的子对象 h=plot(m,sin(m),'r'); set(h,'LineStyle','.','Marker','*','MarkerSize',20); pause(0.001); end 示例 5：普通三维坐标下一个点沿着曲线移动的轨迹。 解： %三维螺旋线坐标 t1=10*pi*(0:;32 x1=cos(t1); y1=sin(t1); z1=-t1; %绘制曲线下面的水平直线 t2=(0:10)/10; x2=x1(end)*(1-t2); y2=y1(end)*(1-t2); z2=z1(end)*ones(size(x2)); %绘制垂直直线 t3=t2; z3=(1-t3)*z1(end); x3=zeros(size(z3)); y3=x3; %绘制曲线上面的水平直线 t4=t2; x4=t4; y4=zeros(size(x4)); z4=y4; %绘制曲线 x=[x1,x2,x3,x4]; y=[y1,y2,y3,y4]; z=[z1,z2,z3,z4]; plot3(x,y,z,'b','linewidth',3) axis off %绘制移动的点 h=line('Color',[1 0 0],'Marker','.','MarkerSize',40,'EraseMode','xor'); n=length(x); i=1;j=1; while 1 set(h,'xdata',x(i),'ydata',y(i),'zdata',z(i)); drawnow %使 Matlab 暂停目前的任务序列而去刷新屏幕 pause(0.0005); f=getframe(gcf); i=i+1; j=j+1; if i&n & j&450 end end33 5、Matlab 在高等数学中的应用5.0 Matlab 中的符号表示1、建立符号变量 Matlab 提供了 sym 和 syms 两个建立符号对象的函数 1）sym 函数，用来建立单个符号变量，调用格式为： 符号变量名=sym(?符号字符串?) 注：符号字符串可以是 常量、变量、函数或者表达式。 2）syms 函数，依次可以定义多个符号变量，调用格式为： syms 符号变量 1 符号变量 2 …… 符号变量 n 注意 符号变量名之间不能使用任何标点符号，只能用空格隔开。 例如： f=sym(?cos(x)?),f=sym(?sin(x)^2=0?); f=sin(x)+cos(x), 2、建立符号表达式 ? ? ? 利用单引号来生成表达式 利用 sym 函数建立符号表达式 使用已经定义过的符号变量组成符号表达式 注意：符号表达式包括符号函数和符号方程，区别在于是否带有等号。 3、符号表达式的运算 (1)四则运算 四则运算分别使用 symadd, symsub,symmul 和 symdiv 来实现加减乘除运算，使用 sympow 实现幂运算。 （注：6.5 以后的版本中已经没有了这些函数，而直接使用+ - * / 即可） (2)提取分子和分母 如果符号表达式是一个有理分式或可以展开为有理分式， 可利用 numden 函数来提 取符号表达式中的分子或者分母，调用格式为：34 [n,d]=numden(s) [n,d] = numden(sym(4/5)) returns n = 4 and d = 5. [n,d] = numden(x/y + y/x) returns (3)因式分解与展开 1）factors(s)分解因式 2）expand(s)展开 3）collect(s)合并同类项 4）collect(s,v)按变量 v 进行合并同类项 (4)符号表达式的化简 1）simplify(S)应用函数规则对 S 进行化简 2）simple(S) 进行综合化简 (5) 符号表达式与数值表达式之间的转换 1）利用函数 sym 可以将数值表达式表示成符号表达式 2）numeric 或者 eval 函数可以将符号表达式转换成数值表达式 3）函数 digits 和 vpa 配合替换函数 subs 进行转换。 digits 函数，digits(D)函数设置有效数字个数为 D 的近似解精度。 vpa 函数 ，R=vpa(S)符号表达式 S 在 digits 函数设置下的精度的数值解。vpa(S,D) 符号表达式 S 在 digits(D)精度下的数值解。 subs 函数，subs(S, NEW, OLD) 4、符号函数的相关运算 （1）复合函数运算。 compose 函数 （2）反函数运算。 finverse 函数 5、符号代数方程求解 (1)线性方程组的求解，函数 linsolve，solve，可以得到方程的精确解 (2)多变量的非线性方程的符号解法，使用函数 fsolve，调用格式有： ? ? ? ? ? X=fsolve(?fun?,X0) X=fsolve(?fun?,X0,options)options 为选择参数输入向量 X=fsolve(?fun?,X0,options,?gradfun?),gradfun 为输入函数在 X 处的偏导数 X=fsolve(?fun?,X0,options,?gradfun?,P1,P2,…)P1,P2 为问题定性参数 [X,options]=fsolve(?fun?,X0,…)返回使用的优化方法的参数 n = x^2+y^2 , d = y*x35 注意：fun 必须使用@指定 示例：定义一个函数 function F = myfun(x,c) F = [ 2*x(1) - x(2) - exp(c*x(1)) -x(1) + 2*x(2) - exp(c*x(2))]; end 求解该函数： c = -1; %定义输入参数 x = fsolve(@(x) myfun(x,c),[-5;-5]) 6、多项式的表示方法 （1）多项式的表示方法――转化为向量问题 对于多项式P( x) ? a0 x n ? a1 x n ?1 ? ... ? an?1 x ? an用行向量表示： P ? [a , a ,..., a , a ] 0 1 n ?1 n ①系数向量直接输入法， MATLAB 自动将向量元素按降幂顺序分配给各系数值。 函数 poly2sym 可以将向量表示的多项式转化为符号多项式表示。 注：由特征多项式生成 的多项式首项系数一定为 1；n 阶矩阵一般产生 n 次多项式。 ③由根创建多项式，由函数 poly 实现。注：若要生成实系数多项式，则根中的复数必 定对应共轭；生成的多项式向量包含很小的虚部时可用 real 命令将其过滤掉。 （2）多项式的运算 ①多项式求值。 输入变量值代入多项式计算时以数组为单元的使用函数 polyval（对应 元素计算）;以矩阵(必须为方阵)为计算单元求多项式的值用函数 ②多项式求根。两种方法，一种是调用函数 roots，另一种是通过建立多项式的伴随矩 阵再求其特征值的方法得到多项式的所有根。 （使用 compan 和 eig 函数） ③多项式的乘除法运算。乘法使用函数 conv（向量卷积） ，除法使用函数 deconv ④多项式微分。微分函数 polyder ⑤多项式拟合。两种方法，一种是由矩阵的除法求解超定方程来进行，另一种是用拟 合函数 polyfit,调用方式为 polyfit(X,Y,n)和[p,s]=polyfit(X,Y,n)36 5.1 求极限Matlab 求极限的命令为 limit(expr, x, a) limit(expr, a) limit(expr) limit(expr, x, a, 'left') limit(expr, x, a, 'right') 其中 limit(expr, x, a)表示求符号表达式 expr 关于符号变量 x 趋近于 a 时的极限， limit(expr)求表示缺省变量趋近于 0 时的极限。 例 15 求下列表达式的极限解 （1）clc, clear syms x f=(3*sin(x)+x^2*cos(1/x))/((1+cos(x))*log(1+x)); s=limit(f,x,0) s=simplify(s)（2）clc, clear syms x a s1=limit(((x+2*a)/(x-a))^x,x,inf) s2=limit(((x+2*a)/(x-a))^x,x,-inf)5.2 求导数Matlab 的求导数命令为 diff(expr) diff(expr, v) diff(expr, sym('v')) diff(expr, n) diff(expr, v, n)37 diff(expr, n, v) 其中 diff(expr)表示求表达式 expr 关于默认变量的 1 阶导数，diff(expr, v, n)和 diff(expr, n, v)都表示求表达式，expr 关于符号变量 v 的 n 阶导数。解 （1）clc, clear syms x y=log(x+sqrt(1+x^2)); s=diff(y,x,2) s=simple(s) %对符号函数进行化简 （2）a=[1,0.5,3.5,6]; da=diff(a)5.3 求一元函数极值计算的 Matlab 程序如下 syms x y=x^3+6*x^2+8*x-1; dy=diff(y); dy_zero=solve(dy), dy_zero_num=double(dy_zero) %变成数值类型 ezplot(y) %符号函数画图5.4 求积分1. 求不定积分 Matlab 求符号函数不定积分的命令为 int(expr) int(expr, v)38 解 syms x I=int(1/(1+sqrt(1-x^2))); pretty(I) %写成符号表达式的形式 2. 求定积分 （1）求定积分的符号解 Matlab 求符号函数的定积分命令为 int(expr, a, b) int(expr, v, a, b)解 syms x I=int(cos(x)*cos(2*x),-pi/2,pi/2) （2）求定积分的数值解 例 20 求下列积分的数值解输入 I=quadl(@(t) (t-3*t.^2+2*t.^3).^(-1/3),eps,0.5) 得到 I=1.4396。 注：quad()函数使用 simpson 法，即用二次曲线逼近被积函数，quadl()采用更高阶 的逼近方法。它们的调用方法是：q=quad(fun,a,b,tol)或者 quadl(fun,a,b,tol)，fun 可以是39 符号方程，也可以是匿名函数。 ）输入 I=dblquad(@(x,y)sqrt(1-x.^2-y.^2).*(x.^2+y.^2&=x),0,1,-0.5,0.5) 得到 I= 0.6028。 注： Matlab 中计算二重积分的函数形式为 q=dblquad(fun,xmin,xmax,ymin,ymax,tol)， 一般只用于积分区域为矩形的情况。为了使用该函数，有时需要将非矩形区域化成矩形 区域。2 2 示例：计算二重积分 ?? ?( x ? y )dxdy 积分区域 ? 由 x=1,y=x 及 y=0 围成。 ??解：首先绘制积分区域，Matlab 程序如下： clear,clc fill([0,1,1,0],[0,0,1,0],'y');hold on%绘制积分区域 fill([0.55,0.6,0.6,0.55,0.55],[0,0,0.6,0.55,0],'r') %绘制单元条 gtext('y=x'); gtext('x=1');pause gtext('y=0') 按照矩形区域调用 dblquad 函数，程序如下： I=dblquad(?(x.^2+y.^2).*(y-x&0)?,0,1,0,1) 求得 I=0.3333 （3）输入 fun3=@(x,y,z)z.^2*log(x.^2+y.^2+z.^2+1)./(x.^2+y.^2+z.^2+1).*(z&=0 &... z&=sqrt(1-x.^2-y.^2)); %该语句使用了续行符... I=triplequad(fun3,-1,1,-1,1,0,1) 得到 I= 0.1273。 %三重积分计算函数40 Matlab 符号求解的程序如下 clc, clear syms r z theta I1=int((r^2+z^2)*r, z, r^2,sqrt(2-r^2)); %求最内层积分 I2=int(I1,r,0,1); %求中间层的积分 I=int(I2,theta,0,2*pi) %求最终的积分结果 pretty(I) %分数线中间显示的格式 Inum=double(I) %把符号解化成数值解 解法二 求数值解的 Matlab 程序如下 clc, clear f=@(x,y,z) (x+y+z).^2.*(z&=x.^2+y.^2 & x.^2+y.^2+z.^2&=2); I=triplequad(f,-sqrt(2),sqrt(2),-sqrt(2),sqrt(2),0,sqrt(2)) 求得积分的值为 2.4486。 注：quad 与 int 的区别，int 可以求具有解析解的不定积分和定积分，quad 只能求解 定积分。5.5 求解微分方程（组）1、微分方程(组)的解析解41 在 MATLAB 中用 D 表示导数， 例如 Dy 表示 y?， D2y 表示 y??， Dy(0)=5 表示 y?(0)=5 等，符号常微分方程求解通过函数 dsolve 实现，调用格式为： dsolve(?e?,?c?,?v?)包括微分方程，初始条件和指定变量三个部分，求解常微分方程 e 在初始条件 c 下的特解，v 描述方程中的自变量，省略则默认以 t 为自变量，若没有初 始条件，则获得通解。 例 1：求 du ? 1 ? u 2 的通解 dt解：Matlab 程序为： dsolve('Du=1+u^2','t') ?d2 y dy ? 4 ? 29 y ? 0 ? 的特解 例 2：求微分方程 ??dx 2 dx ? ??y (0) ? 0, y?(0) ? 15 解：Matlab 程序为： dsolve('D2y+4*Dy+29*y=0','y(0)=0,Dy(0)=15','x') ? dx ? dt ? 2 x ? 3 y ? 3 z ?? ? dy 例 3：求常微分方程组的通解 ? ? 4 x ? 5 y ? 3 z ??dt ? dz ? ? 4x ? 4 y ? 2z ??dt 解：Matlab 程序为： [x,y,z]=dsolve('Dx=2*x-3*y+3*z','Dy=4*x-5*y+3*z','Dz=4*x-4*y+2*z') 一阶齐次线性常微分方程组和非齐次线性方程组可以直接通过矩阵操作求解。 对于一阶齐次线性微分方程组： X ? ? AX , X ? ? x1 ,xn ? , A ? (aij )n?n 其解的形式为：X (t ) ? e A(t ?t0 ) X 0对于非齐次线性方程组： X ? ? AX? f (t ), X (t 0 ) ? X 0 解的形式为：X (t ) ? e A(t ?t0 ) X 0 ???? e A(t ? s ) f (s)ds? ? ? ?? 例 4： X ? ? ? 0 2 ?1? X , X (0) ? ?2 ?? ? ?0 0 3 ? ?1 ?2 1 3? ?1 ??42 解：a=[2,1,3;0,2,-1;0,0,2]; x0=[1;2;1];x=expm(a*t)*x0 pretty(x); %化简结果?1 0 0 ? ? 0 ? ? 0 ?? ? ? ? ? ? ?? 例 5： X ? ? ? 2 1 ?2 ? X ? ? 0 ? , X (0) ? ?1 ?? t ? ? ? ?3 2 1 ?e cos(t ) ?1解： a=[1 0 0;2 1 -2;3 2 1]; fs=[0;0;exp(s)*cos(2*s)]; x0=[0;1;1]; tx=int(expm(a*(t-s))*fs,s); %先求不定积分，再计算定积分的结果，提高速度 xstar=subs(tx,s,t)-subs(tx,s,0); x=expm(a*t)*x0+ x=simple(s),pretty(x) %simple 直接显示化简结果，simplify 显示化简过程，pretty 修 改显示格式 2、微分方程（组）的数值解 当难以求得微分方程的解析解时，可以求其数值 解。Matlab 中求解微分方程数值解 的方法有以下几个：ode45，ode23，ode113，ode23s，ode15s。调用形式如下： [t,x]=solver(‘f’, ts, x0, options)其中，t 是自变量的值，x 是函数值，solver 表示上述 5 种求解器的一种；‘f’是由待解方程写成的 m 文件名，ts=[t0, tf]表示自变量的初值和 终值，options 用于设定误差限，通过函数 options=odeset(‘reltol’, rt, ‘abstol’, at)分别设 定相对误差和绝对误差。 不同的函数代表不同的内部算法，ode45 表示 4/5 阶龙格-库塔-费尔贝格算法（RK 算法），是解非刚性常微分方程的首选方法，ode23 采用 2/3 阶 RK 方法，ode113 采用 多步法，效率一般比 ode45 高。ode15s，ode23s,ode23t,ode23tb 用来解刚性常微分方程。 注 1：在解 n 个未知函数的的方程组时，x0 和 xn 均为 n 为向量，m 文件中的待解 方程组应以 x 分量形式写出； 注 2：使用 Matlab 求数值解时，高阶微分方程必须等价的变换成一阶微分方程组。 ? ?d2 x 2 dx ? ? ? 1000(1 ??x ) x 0 ? 例 1：求解 van der Pol 刚性微分方程 ??dt 2 的数值解 dt ? ??x(0) ? 2; x '(0) ? 043 解：令 y1=x，y2=y1?，则微分方程变为一阶微分方程组： ? y1' ? y 2 ? 2 ??y 2 ' ? 1000(1 ? y1 ) y 2 ? y1 ? y1(0) ? 2; y 2(0) ? 0 ?? 写成 m 文件如下所示： function dy=fun1(t,y) dy=zeros(2,1); dy(1)=y(2); dy(2)=1000*(1-y(1)^2)*y(2)-y(1); 取 [t0 ,tf]=[0,3000] ，输入命令： [t,y]=ode15s('fun1',[0,3000],[2,0]); plot(t,y(:,1),'-') ? y1 ' ? y2 y3 ?? ? ? y1 y3 ? y2 ' 例 2：求解微分方程组 ?? ? y3 ' ? ?0.51 y1 y2 y1 (0) ? 0, y2 (0) ? 1; y3 (0) ? 1 解：建立方程组的 m 文件： function dy=fun2(t,y) dy=zeros(3,1); %初始化 dy(1)=y(2)*y(3); dy(2)=y(1)*y(3); dy(3)=-0.51*y(1)*y(2); 取 t0=0，tf=12，Matlab 命令如下： [t,y]=ode45('fun2',[0 12],[0 1 1]) plot(t,y(:,1),'-',t,y(:,2),'*',t,y(:,3),'+')?? x ? ? ( y ??x) ?? y ? ? x ? y ? xz 在初始条件[x0,y0,z0]=[5,13,17]时的解，其 例 3：求解 Lorez 方程 ?? ?? z ? ?? z ??xy ??中 ? ? 10, ? ? 28, ? ? 8 / 3 。 解：rho=10;beta=29;lamda=8/3;44 f=@(t,Y)[rho*(Y(2)-Y(1));beta*Y(1)-Y(2)-Y(1)*Y(3);-lamda*Y(3)+Y(1)*Y(2)];%方程 组必须使用列向量表示 [t,y]=ode45(f,[0,30],[5,13,17]) subplot(2,2,1) plot(t,y(:,1),'*') %x 曲线 subplot(2,2,2) plot(t,y(:,2),'X') subplot(2,2,3) plot(t,y(:,3),'O') subplot(2,2,4) plot3(y(:,1),y(:,2),y(:,3)) %绘制空间轨迹图 3、边值问题的 Matlab 数值解 Matlab 中使用 bvp4c（bvp5c）和 bvpinit 命令求解常微分方程中的两点边值问题， 其调用格式为： sol = bvp4c(odefun,bcfun,solinit,options) solinit = bvpinit(x,yinit,parameters) %该函数用来给 sol 提供初始猜测解。5.6 级数求和Matlab 级数求和的命令为 r = symsum(expr, v) r = symsum(expr, v, a, b) 其中 expr 为级数的通项表达式，v 是求和变量，a 和 b 分别为求和变量的起始点和 终止点，若没有指明 a 和 b，a 的默认值为 0，b 的默认值为 v-1。 无穷数列的累加称为级数，当取其前面若干有限项时，得到的是部分和，数列 a 累 加形成的新序列可用 s=cumsum(a)实现。 a 的长度为 n， 则 s 的长度也为 n， 即每一个 s(k) 是数组 a 中前 k 项的和。注意 cumsum 与 sum 的区别，sum(a)得到的是一个数，即序列 s 中最后一项。 示例： clc,clear k=1:10; a=1./k.^2;45 s=cumsum(a) ss=sum(a) 例 22 求如下级数的和解 （1）syms n f1=(2*n-1)/2^n; s1=symsum(f1,n,1,inf) （2）syms n f2=1/n^2; s2=symsum(f2,n,1,inf) 补充一点： Matlab 求数组元素乘积的函数 B=prod(A) 如果 A 是一个 m 行一列的 （向 量），则这种用法即返回这 m 个元素的乘积； 如果 A 是一个 m 行 n 列的矩阵，则 A 的每一列都被看做一个 m 行 1 列的向量，分别计算每个向量中元素的乘积，返回给 B， 因此 B 是一个 1 行 n 列的数组。5.7 Matlab 求解优化问题Matlab 求解各种形式的极值问题被集成在优化工具箱中。 下面是优化工具箱 4.0 版本的功能示意图：46 47 变量 f描述 线性规划目标函数 f*X 或二次规划的目标 函数 X’*H*X+F*X 中线性项的系数向量 非线性优化的目标函数. fun 必须为行命令 对象或 M 文件、 嵌入函数、 或 MEX 文件的 名称 二次规划目标函数 X’*H*X+f*X 中的二次 项的系数矩阵 A 矩阵和 b 向量分别为线性不等式约束： AX≤b 中的系数矩阵和右端向量 Aeq 矩阵和 beq 向量分别为线性等式约束 Aeq*X=beq 中的系数矩阵和右端向量 X 的下限和上限向量 迭代初始点坐标 函数最小化的区间 优化选项参数结构调用函数 linprog,quadprog fminbnd,fminsearch,fminunc, fmincon,lsqcurvefit,lsqnonlin, fgoalattain,fminimax quadprog linprog,quadprog,fgoalattain, fmincon,fminimax linprog,quadprog,fgoalattain, fmincon,fminimax linprog,quadprog,fgoalattain, fmincon,fminimax,lsqcurvefit, lsqnonlin 除 fminbnd 外所有函数 fminbnd 所有优化函数funHA,bAeq,beqvlb,vub X0 x1,x2 options变量描述 描述退出条件： exitflag&0 表示目标函数收敛于 x exitflag=0 表示已达到函数评价或迭代的最大次 数 exitflag&0 表示目标函数不收敛 由优化函数求的 x 值调用函数exitflagx所有优化函数 linprog,quadprog,fgoalattain fmincon,fminimax,lsqcurvefit lsqnonlin,fminbndfeval解 x 处的目标函数值output包含优化结果信息的输出结构 所有优化函数 iterations,Algorithm,FuncCount（函数评价次数）进行极值计算时有两种方式，通过命令 optimtool 打开优化工具箱 GUI 求解和使用48 上述提供的命令函数进行求解。 示例一：使用单纯形计算最小值 min f ? ?4 x1 ? x2 ? ? x1 ? 2 x2 ? 4 ? 2 x ? 3 x ? 12 ?? 2 s.t. ? 1 ? x1 ? x2 ? 3 x1 , x2 ? 0 方法一：通过优化工具箱求解 方法二：输入 Matlab 命令： clc,clear f=[-4;-1]; A=[-1 2;2 3;1 -1]; b=[4;12;3]; [x,fval,exitflag,output]=linprog(f,A,b) 示例二：非线性优化 f ( x) ? ? x1 x2 x3 ? x ? 2 x2 ? 2 x3 ? 0 ?? s.t. ? 1 ? x1 ? 2 x2 ? 2 x3 ? 72 方法一：通过优化工具箱求解 方法二：输入 Matlab 命令： clc,clear f=@(x)(-x(1)*x(2)*x(3)); A=[-1 -2 -2;1 2 2]; b=[0;72]; [x,fval,exitflag,output]=fmincon(f,[10;10;10],A,b) 其他求解极值的函数主要 有： (1)x=fminbnd(fun,x1,x2)，求函数在区间[x1,x2]上的极小值对应的自 变量值 示例三：求 f(x)=x^4-x^2+x-1 在区间[-2,1]上的极小值。 解： Matlab 命令如下： [x,feval,exitflag,output]=fminbnd('x^4-x^2+x-1',-2,1) (2) x=fminsearch(fun,x0)无约束非线性极小值求解，x0 为起始搜索点。49 示例四：fminsearch 函数求 f(x)=1/((x1-2)^2+3)-1/(2(x2+1)^2-5),初始点[0,0]， (3) x=fminunc(fun,x0,options)求解函数的一个局部极小值 (4) x=fminimax(fun,x0,A,b,Aeq,beq,lb,ub,nonlcon)nonlcon 包括非线性约束; (5) x=lsqnonlin(fun,x0,lb,ub);非线性最小二乘优化 (6) x=linprog(f,A,b,Aeq,beq,lb,ub,x0)求解线性规划问题 (7) x=bintprog(f,A,b,Aeq,beq,x0)求解 0-1 规划 (8) x=quadprog(H,f,A,b,Aeq,beq,lb,ub,x0);求解二次规划5.8 Matlab 求解插值和拟合问题5.8.1 Matlab 插值工具箱 1、一维插值函数 interp1，命令格式为： yi=interp1(x,Y,xi,?method?)其中 method 指定的方法为 nearest,linear,spline,cubic， pchip （分段三次 Hermite 插值） ， 所有的插值方法都要求 x 是单调的。 x 与 Y 是具有相同大小 的向量，求在 xi 点处的插值函数值 yi 一维快速插值函数 interp1q，命令格式为： yi = interp1q(x,Y,xi)，其中 x 是单调递增的列向量，Y 是列向量或者行数为 length(x) 的矩阵，xi 是一个列向量。 yi = interp1(x,Y,xi,method,'extrap') %对于超出 x 范围的 xi 中的分量将执行特殊的外 插值法 extrap。 示例 1：在一 天 24 小时内，从零点开始每间隔 2 小时测得的环境温度数据分别为 12，9，9，10，18 ，24，28，27，25，20，18，15，13，推测中午 12 点（即 13 点）时的温 度． 解： x=0:2:24; y=[12 9 9 10 18 24 28 27 25 20 18 15 13]; a=13; y1=interp1(x,y,a,'spline') xi=0:1/3600:24; %绘制一天内的温度变化曲线 yi=interp1(x,y,xi, 'spline'); plot(x,y,'o' ,xi,yi) 示例 2： x=0:0.12:1; y=(x.^2-3*x+5).*exp(-5*x).*sin(x); plot(x,y,'ro',x,y);50 x1=0:0.02:1; %要插值点 y1=(x1.^2-3*x1+5).*exp(-5*x1).*sin(x1); y2=interp1(x,y,x1);y3=interp1(x,y,x1,'spline'); y4=interp1(x,y,x1,'nearest');y5=interp1(x,y,x1,'cubic'); plot(x1,[y2' y3' y4' y5'],':',x,y,'ro',x1,y1); legend('linear','spline','nearest','cubic','样本点','原函数'); %计算各插值的最大误差 [max(abs(y1-y2)),max(abs(y1-y3)),max(abs(y1-y4)),max(abs(y1-y5))]; 2、二维数据内插值 interp2，插值函数为二元函数。 ZI = interp2(X,Y,Z,XI,YI,?method?) X 和 Y 分别是 m 维和 n 维的向量，表示节点，Z 为 n*m 的矩阵，表示节点值，XI， YI 是为一维数组，表示插值点，XI 与 YI 是方向不同的向量，一个是行向量，一个是列 向量，ZI 是行数为 YI 的维数，列数为 XI 的维数的矩阵，表示得到的插值。若 XI 与 YI 中有在 X 与 Y 范围之外的点，则相应地返回 nan（Not a Number） 。method 为插值 方法。 如果在某区域测量了若干个节点的高程（节点值） ，为了画出比较精确的等高线图， 需要先插入更多的点（插值点） ，然后计算这些点的高程。 如果是三次样条插值，可以使用命令 pp=csape({x0,y0},z0,conds,valconds);z=fnval(pp,{x,y}) 示例 1： years = 0; service = 10:10:30; wage = [150.697 199.592 187.625 179.323 195.072 250.287 203.212 179.092 322.767 226.505 153.706 426.730 249.633 120.281 598.243]; w = interp2(service,years,wage,15,1975) %计算在 1975 年工作了 15 年的员工的工资[x,y]=meshgrid(-3:0.6:3,-2:0.4:2); z=(x.^2-2*x).*exp(-x.^2-y.^2-x.*y); surf(x,y,z);figure [x1,y1]=meshgrid(-3:0.2:3,-2:0.2:2); z1=interp2(x,y,z,x1,y1); z2=interp2(x,y,z,x1,y1,'cubic'); z3=interp2(x,y,z,x1,y1,'spline'); surf(x1,y1,z1);figure %比较误差 z0=(x1.^2-2*x1).*exp(-x1.^2-y1.^2-x1.*y1);51 surf(x1,y1,abs(z0-z2));axis([-3,3,-2,2,0,0.08])surf(x1,y1,abs(z0-z3)); axis([-3,3,-2,2,0,0.08]) interp2 函数能够较好的进行二维插值运算，但是它只能处理以网格形式给出的数 据，如果数据是随机给出的，此时可以使用 griddata 函数。其调用格式为： 格式:z=griddata(x0,y0,z0,x,y,?method?) ? v4 ?:MATLAB4.0 提供的插值算法，公认效果较好； ?linear ?：双线性插值算法（缺省算法） ； ?nearest?：最临近插值； ?cubic?：双三次插值。clear x=-3+6*rand(199,1);y=-2+4*rand(199,1); z=(x.^2-2*x).*exp(-x.^2-y.^2-x.*y); plot(x,y,'*');plot3(x,y,z,'*'); [x1,y1]=meshgrid(-3:0.2:3,-2:0.2:2); z1=griddata(x,y,z,x1,y1,'cubic'); surf(x1,y1,z1); z2=griddata(x,y,z,x1,y1,'v4');surf(x1,y1,z2); %误差比较 z0=(x1.^2-2*x1).*exp(-x1.^2-y1.^2-x1.*y1); surf(x1,y1,abs(z0-z1));axis([-3,3,-2,2,0,0.15]) surf(x1,y1,abs(z0-z2));axis([-3,3,-2,2,0,0.15]) 3、 三维插值运算函数 interp3,n 维网格插值 interpn， 其调用格式同 interp1 和 interp2， 对应的三维网格生成函数为[x,y,z]=meshgrid(x1}
我要回帖
更多推荐
- ·郑莉君社会兼职网官网入口
- ·备考圈高考志愿填报机构哪家好怎么样？
- ·道路与桥梁工程技术专业介绍CAD内容简介
- ·台湾有几个建行银行大楼建筑简介
- ·政府审计的四种类型内容有哪些
- ·蒸汽卧式往复泵泵比较好的去哪买？
- ·未来含金量最高的10大资格证含金量排名有哪些
- ·ListenAndServeTLS函数不能android 连接httpss怎么回事
- ·谁能帮我写一篇作文400字600字的作文？和亲人有关的，最好是父母或者爷爷奶奶，像 回忆我的母亲 那种的
- ·有在使用蓝鹰教育环保的净化设备的吗？感觉如何？
- ·澳大利亚新西兰地图、新西兰如何防止市领导利用职务上的便利为他人谋取利益并收受巨额财物?
- ·专家解读人与技术之间复句到底是什么关系系
- ·objective c中的数组-c以json方式向后台传数组，应该怎么办
- ·我有个形容朋友关系好的成语叫琳珊能有什么成语把俩个字
- ·知道其他投标人资格证明的资格证明文件有误应该怎么质疑
- ·栞遍这个字的组词是什么怎么读，栞怎么念，栞怎么拼音，栞怎么组词
- ·花呗逾期一个月3个月是不是永久不能使用？
- ·简便画桂林森美装饰怎么样拉怎么画
- ·有谁得到过感情老师的帮助挽回感情的短信了感情
- ·matlab求matlab 解微分方程程my"+cy'+ky=F(t)，m,c,k均为常数，当F(t)以离散的表格形式给
- ·整理卧室宪法的感想300字
- ·输出的功率来源于直流电源功率，可不可以是交流电源
- ·06款赛欧06款别克君威水温问题表不准怎么回事
- ·现代领动殡仪车副驾像乔任梁哪里有电流声怎么回事
- ·若5属于，excel求x值的值
- ·雨露计划，教育补助申请，大专学校应该填写高职和大专吗？
- ·练扎桩摇的时候，练钢琴小拇指关节酸痛撑得很疼，是不是姿势错了，求
- ·怀孕才一个多月，就不小心如何让同事认识你知道，会不会对
- ·药流没流干净不想清宫不完全清宫后5天有褐色分泌物是为什么
- ·男孩7岁，时不时鳄鱼张嘴原因深呼吸一下是什么原因
- ·我有乙肝，会抽烟影响试管成功率吗我做试管吗
- ·杨西宁治疗男性狂躁症的治疗方法怎么样
- ·我经常睾丸疼痛会有疼痛腰酸和遗精请问应该用什么药
- ·为什么胰岛素抑制胰高血糖素能口服吗，反过来不行
- ·身上起龟头红疙瘩怎么回事痒痒还发烧是怎么回事
