数值计算与MATLAB教材

数值计算与MATLAB0前言实际问题→数学模型→数值计算方法→程序设计→上机计算→求出结果

“数值计算”、“数值分析”为同义词，这门学科有如下特点：面向计算机：可实现性有可靠的理论基础：理论明确要有好的计算复杂性：可行性可实现数值验证：可验证性算法误差明确、并可满足要求：可信性

主要内容：MATLAB基础（方程组、最小二乘法、回归分析）、MATLAB数据可视化、常用函数、复数、向量计算、插值法、曲线拟和、数值微积分、多项式、级数展开、常微分方程求解、符号计算等，其中有些问题包含在各章节中，主要以MATLAB为线索，进行展开。1MATLAB基础

1.1基本数据类型与运算MATLAB数据类型及常数表达名称：字母开头,大小写敏感，缺省变量：ans

实数：科学计数法，特殊数表达eps:浮点数的最小单位；Inf无穷大；NaN不是一个数值；pi;realmax;realmin复数：a+bi,或a+bj

，a,b为实数,i、j虚数单位矩阵：[………],其中各项之间用空格、逗号、分号分开，分号换行1.1基本数据类型与运算特殊矩阵：eyeonesrandrandnzeros…(n,m)/(n):nbym,nbyn线性向量：线性a:k:b

产生向量示例：1）用360个弧度表示一圆周；

2）产生200*200个元素为（0，1）正太分布随机数矩阵；1.2运算参与运算的基本单位：矩阵单个实数、复数、向量都是矩阵，都可以参加运算。矩阵运算：+-*/\^数组运算：

.*./.^其中矩阵的除法

/:B/A——B*A-1\:A\B——A-1B

如果XA=B则:X=B/A

如果AX=B则:X=A\B‘转置.’非共轭转置矩阵与常数运算示例生成200*200个元素为复数元素矩阵，复数实部与虚部均为（10，5）正态分布随机数产生2000个角度为一周均匀分布，直径为正态分布的点云，方差为10，圆心为（0，0）点。exp求解1.3矩阵基本操作矩阵的转置矩阵元素访问及矩阵重构A（向量，向量）几个常用函数

fliplr

flipudrot90repmat(A,m,n)

reshape(A,m,n)示例：

A:100*100,B:40*40,将B嵌入到A中心位置

A:300*400,B:300*400,合成“块交错”的C，块的大小为5*5，奇数块来至A,偶数块来至于B1.4线性方程组、最小二乘法与回归分析线性方程一般形式AX=BX=A\BA、B可以是复数？复数方程，可以做！A列数与X行数，B行数一致，未知数与方程数A行数少于列数A行数多于列数线性方程一般形式AX=B1）B列数大于1，相当于一个A对应多个B的方程2）A行数小于列数，相当于方程数小于未知数个数3）A行数大于列数，相当于方程数多于未知数个数1.4线性方程组、最小二乘法与回归分析未知数多，方程数少：得到一个基本解未知数少，方程数多：采用最小二乘法，得到一个最综合的解，也就是使得“=”左端计算的值，与B的差最小。最小二乘法：对于y=ax+b

形式的线性方程，如果已知点(x1,y1)和(x2,y2)满足方程，则同时满足y1=a*x1+b;y2=a*x2+b，就可以建立方程即可求出(a,b)值,如果多于2个点，则应该找到一条直线与点最近，实际上是求出线到点的y距离的总和最小，数学上利用y距离的平方和确定总误差。即E=SUM(y*-y)2,E中仅包含两个未知数a,b.求E最小可以利用偏导等于0的方程。1.4线性方程组、最小二乘法与回归分析X=1:5’X=[xoese(size(x))]Y=[1030324548]’A=x\y，即为线性方程的两个系数误差多大？怎么求？（相对误差、绝对误差）

如何做非线性回归x=1:10;y=[20212424.62630.129.7384353]方程形式为y=ax2+bx+c实际上1.4线性方程组、最小二乘法与回归分析MATLAB实现:x,y

转换为列向量X=[x.^2xones(size(x))]A=x\y,a(1)->aa(2)->ba(3)->c如何理解ones(size(x))?Size返回什么其他方程形式回归！1.5环境管理Clc：清楚控制台（操作窗口）Clear：清楚内存变量，其他形式？清楚一个（一批）特定变量help内容QuitSAVEfnameXYZ-ASCII-DOUBLeload(filename)whowhos输出格式设定formatshort/long/shorte/longe/shortg/longg/formathex/bank/rat2常用函数与数据可视化2.1常用函数

Magic(n)魔方阵pascal(n)阵hadamard(n)n为2的密正交阵

calendar(yyyy,mm)月历；clock时间向量

log2log10logsqrtexpabsimagrealangleconjfix向0取整floorceilround舍入

sign数单位值

Sincostanasin

acos

atanatan22.2平面曲线绘制plot(y)

如果y为实数向量则绘制各点折线，复数向量以实部为x,虚部为y绘制曲线，y为矩阵则以列为线绘制折线plot(X1,Y1,...)plot(X1,Y1,LineSpec,...)绘制（x1,y1）线并指定线的特性线型：-实线--划线:点线-.点划线数据标志：ox+sd*.^v><ph

线的颜色：ym(agenta)c(yan)rgbw(blac)k

2.2平面曲线绘制绘制120度差的三条sin曲线；绘制多边形；曲柄连杆机构轨迹；绘制动刀轨迹平面几何变换：点阵的旋转、缩放、平移2.2平面曲线绘制Bar(x,y,width)Polar(theta,rho)Pie(x,explode)SemilogxSemilogyLoglog示例：绘制螺旋线；绘制传递函数的幅频特性曲线，相频特性曲线（利用半对数曲线)2.3视图管理与控制1）坐标轴axis([xmin

xmax

ymin

ymax

zmin

zmax])axisautoaxismanualaxistightaxisfillaxisequalaxissquareaxisoffaxison2）分割/设定绘图区

Subplot(m,n,p)分割成m行n列，当前设为p列3)Holdon/off4）titlexlabel

ylabel

zlabel

指定标题平面曲线示例3数据分析k=convhull(x,y):包围所有点(x,y)的凸多边形的提取，返回值为向量in=inpolygon(X,Y,xv,yv)多边形内点的检测,返回0,13.1常规统计分析常规分析f=factor(n)返回n整数的素数因子primes(n)

返回小于等于n的所有素数P=perms(v)

所有排列sum,mean,std,median,mode,:

格式：函数（a,dim）a为矩阵dim=2为列间计算max,min

格式：函数（a,b,dim）a,b相同大小[n,x]=hist(y,m):m为组数，n为个组计数，x为组中值polyarea(X,Y)多边形面积

3.2有限差分Diff(x):差分Cumsum（x):累加.求导、积分示例：一点，在X轴上运动，每3秒观察一次，数据为03255648108912849493，求速度3.3FFT变换傅里叶积分变换实质：FFT(x):返回值包含两部分，幅值与相位，并且是左右对称。示例：s=10*sin(t)+7*sin(3*t+0.8)+5*sin(5*t+1.2)+3*sin(8*t+1.3)3.4数值滤波-IIF滤波infiniteimpulseresponse数值滤波处理y(n)=b(1)*x(n)+b(2)*x(n-1)+...+b(nb+1)*x(n-nb)-a(2)*y(n-1)-...-a(na+1)*y(n-na)y=filter(b,a,X)其中a的第1项一定为1如：对一个序列向前3项平均处理:x=[001111]y=filter([111]/3,1,x);结果为：000.330.66711*指数平滑：y(k)=ax(k)+(1-a)y(k-1)y=filter(0.7,[1-0.3],x);*声音的ECHO处理，设回声时差为0.3s?,回声强度0.6?，采样频率320004多项式与插值4.1多项式表达采用实数行向量表示降幂多项式系数，如：p=[3425]表示为4.2多项式运算与求解

conv(p1,p2)多项式相乘，卷积计算

roots(p)求解多项式的根

polyval(p,x)多项式求值

polyfit(x,y,n)拟和多项式Poly(a):如果a为矩阵则返回矩阵A的特征多项式p；如果a为向量，返回以向量中的元素为根的多项式。[q,r]=deconv(u,v)：u/v，q商多项式，r余数多项式4多项式与插值4.3插值

yi=interp1(x,Y,xi,method)一维插值，method可为：

'nearest‘;'linear';'spline';‘cubic’

*插值：利用插入点X附近的数据，按一定方法求出一个表达式，再利用表达式计算该点

spline:样条插值；cubic:立方插值

zi=interp2(X,Y,Z,XI,YI,metho):二维插值示例：数据插值、曲面插值（10个值）插值示例Sinsin面z=sin(sqrt(x.^2+y.^2)/2)4多项式与插值4.4

数值网格化[XI,YI,ZI]=griddata(x,y,z,xi,yi)

。x,y,z为实测数值，xi,yi为网格坐标，返回网格坐标值，如：x=rand(300,1)*4-2;y=rand(300,1)*4-2;z=x.*exp(-x.^2-y.^2);[xi,yi]=meshgrid(-2:0.2:2);zi=griddata(x,y,z,xi,yi);5三维与特殊形式绘图5.1三维曲线

plot3(X1,Y1,Z1,LineSpec,...)：三维空间曲线

Bar3(x,y,width)

cylinder(r,n)

5.2三维曲面1）产生网格数据[X,Y]=meshgrid(x,y)2）绘制网格曲面

mesh(X,Y,Z)mesh(Z)mesh(...,C)meshc(...)meshz(...)3）绘制曲面surf(Z)surf(X,Y,Z)surf(X,Y,Z,C)surfc(...)4）绘制等值线contour(Z,n)r=sqrt(x.^2+y.^2)+0.1;z=sin(r)./rPeaks(N)5.3曲面平滑与颜色

view(az,el)：观察方向1）平滑shadingflat/faceted/interp

2）利用colormap(map)

设定曲面颜色，map为64by3的矩阵，为颜色索引，系统提供的map为：AutumnbonecoolcopperflagGrayhotjetlinesprismspringSummerwinter3）Colormap([r,g,b]);6逻辑函数逻辑关系：

<><=>===~=运算：&|~All(a)TesttodetermineifallelementsarenonzeroAny(a)Testforanynonzeros

k=find(x);[i,j]=find(X)Findindicesandvaluesofnonzeroelements示例1）将a中大于0的元素加和2）将a中小于0的元素加57程序设计基础7.1函数的定义

Function[mean,stdev]=stat(x)n=length(x);mean=sum(x)/n;

stdev=sqrt(sum((x-mean).^2/n));其中：Function[mean,stdev]=stat(x)为函数首部，“=”左侧为返回值形式变量，右侧为函数名，括号内为入口的形式参数；其中可以用“%”引导注释内容。文件名应和函数名相同。7.2程序结构控制分支结构

ifexpressionstatementselse

statementsend开关结构switch表达式（标量或字符串）

case值1

语句组ACase值2

语句组B…….Otherwise

语句组Nend扫描循环

forvariable=expressionstatementsendExpression：向量逻辑循环whileexpressionstatementsend示例：返回处入口向量升序、及降序、及相邻项相等的个数function[abc]=myf1(x)l=length(x);%%求x中元素数量

a=0;b=0;c=0;fork=1:(l-1)ifx(k)>x(k+1)a=a+1;elseifx(k)<x(k+1)b=b+1;elsec=c+1;endendend常规程序设计方法太复杂，如何运用MATLAB的功能？设计出更简洁程序？function[abc]=myf1(x)

dx=diff(x);a=sum(dx>0);b=sum(dx<0);c=sum(dx==0);EndFuncitony=myf1(x);

dx=sign(diff(x))+2y=accumarray(dx',ones(size(dx)));End;十进制乘法计算，可以使用卷积，如何十进制进位，字符操作？7.3几个函数Accumarray(s,v):s列向量，v行向量，s,v元素数量相同根据s指定，累加v，生成新的向量，其中s元素指出累加的位置，如：v=11:15,s=[21241]’,结果为：27；24；0；14circshift(A,shiftsize)：循环移动矩阵的行和列，shiftsize

：[l,r]l：移动行数，r：移动的列数，正：右下；负：左上

a=1:5;a=a’;b=circshift(a,2),结果为：[45123]’Sort(a,dim)Colormap(v):v是个n*3向量；hot;gray;cool等，表示colorindexImage(a)：以colorindex形式显示矩阵Globalv1v2Clf：clearfigurev=input(‘prompt')Ginput(n)Pause(d):控制程序停止ds，d可以是小数程序设计示例分析显示移动的彩色块移动的SIN函数旋转的多边形曲柄连杆机构的运动一种“绘画”技术functiona=test1;n=input('Interthelengthof1thline');

clf;holdonaxisequal;axisoffsubp(0,n*exp(i*pi/2));endfunctionb=subp(p,l)

plot([p;p+l]);ifabs(l)>5subp(p+l,0.71*l*exp(i*pi/4));subp(p+l,0.71*l*exp(-i*pi/4));End;end8微分方程数值解与数值积分8.1微分方程数值解

可以计算出一个序列的y，即为数值解，可以使用Cumsum(x)函数实现累加过程，但先要给定初值。例如：y’=cos(x),y=?EulerR-K方法微分方程数值解实际上是前后两点斜率的平均值。上式叫二阶Runge—Kutta法上式叫四阶Runge—Kutta法MATLAB微分方程数值解法[t，x]=ode23(f，ts，x0，options)[t，x]=ode45(f，ts，x0，options)分别为二阶，四阶R-K方法，其中f微分方程（组）的描述函数；ts:为微分点的时间序列，求出的x,为与ts对应的解序列，x0为初值，options为求解过程误差控制参数，可以省略示例11）定义方程函数

functiondx=mf(t,x);

dx=cos(t);end2）求解

ts=0:0.01:10;[t,x]=ode45(‘mf’,ts,0);Plot(x);%%显示曲线1）定义方程函数

functiondy=mf(t,y);

dy=-t*y.^2;end2）求解

ts=0:0.01:1;[t,y]=ode45(‘mf’,ts,2);Plot(y);%%显示曲线示例2其中：a=20;b=40;c=15;t的范围[0,0.5],x(0)=0;y(0)=0注意：为二元微分方程1）编辑函数functiondx=mf(t,x)a=20;b=40;c=15;s=sqrt((c-x(1))^2+(a*t-x(2))^2);

dx=[b*(c-x(1))/s;b*(a*t-x(2))/s];end函数返回值[dx/dt;dy/dt]2）ts=0:0.01:0.5;x0=[00];

[t,x]=ode45(‘mf’,ts,x0);%%采用缺省的精度求解X(1)——x,x(2)——y示例31)定义方程函数functiondu=mf(t,u);du=[-10*u(1)+9*u(2);10*u(1)-11*u(2)];end2)求解

ts=0:0.01:0.5;[t,u]=ode45(‘mf’,ts,[12]);示例4定义微分方程函数functiondy=mf(t,y);

dy=[0;0];dy(1)=y(2);dy(2)=(7*y(2)-5*y(1)+10)/2;end求解：[t,y]=ode45('mf',0:0.01:1,[21])实际上y(2)’，应该与y(1)相同微分simulink解法上个方程8.2数值积分trapz(x,y)：梯形法求积分，求面积q=quad(@fun,a,b)：

Simpson法积分示例：定义被积函数functionyy=mf(x);

yy=1./(x.^3-2*x-5);end积分运算Q=quad(@mf,0,2)9函数极值与优化问题fplot(fun,limits):绘制函数曲线如：fplot(@sin,[0,2*pi]);x=fminbnd(fun,x1,x2):返回单变量函数fun，在[x1,x2]区间上的最小值，如：fminbed(@sin,-1,1),结果为-1x=fminsearch(fun,x0)：返回多变量函数fun，在x0附近的最小值如：functionf=myf(x);f=100*(x(2)-x(1).^2).^2+(1-x(1)).^2;endx=fminsearch(@myf,[01])9函数极值与优化问题x=bintprog(f,A,b,Aeq,beq)：求多变量函数线性函数f，(1/0)规划的极小值，约束条件：AX<=b;AeqX=beq;X取值范围1/0,其中f为向量表示的线性目标函数如：求f(x)=-9x1-5x2-6x3-4x4极小值，满足

6x1+3x2+5x3+2x4<=9;x3+x4<=1;*x=fgoalattain(fun,x0,goal,weight,A,b,Aeq,beq,lb,ub)：目标规划x=fmincon(fun,x0,A,b,Aeq,beq,lb,ub)：条件极值9函数极值与优化问题x=fsolve(fun,x0,options)：返回fun在x0的解，如：在[-5-5]附近的解，定义函数

functionf=myf1(x)f=[2*x(1)-x(2)-exp(-x(1));-x(1)+2*x(2)-exp(-x(2))];end求解：x=fsolve(@myf1,[-5;-5])options=optimset('param1',value1,'param2',value2,...)Display:Levelofdisplay.'off'displaysnooutput;'iter'displaysoutputateachiteration9函数极值与优化问题Fsolve示例分析：四杆机构求解已知量：l1,l2,l3,l4;q1,q4未知量：q2,q3function[bc]=jg();globall1l2l3l4q1q4;l1=20;l2=40;l3=30;l4=50;q4=0;x0=[11];b=[];c=[];op=optimset('Display','off');forq1=0:pi/180:2*pi;x=fsolve(@my,x0,op);x0=x;b=[b;l1*exp(i*q1)];c=[c;l2*exp(i*x(1))];endendfunctionf=my(x)globall1l2l3l4q1q4;%%x(1)-q2;x(2)-q3f=[l1*cos(q1)+l2*cos(x(1))+l3*cos(x(2))-l4*cos(q4);l1*sin(q1)+l2*s