昆明理工大学上机安排2-数值计算

内容：本次上机主要练习常见的数值计算问题，包括线性代数、函数分析、数值微积分、常微分方程等，重点练习如何利用MATLAB提供的函数来实现数值计算，对于数学理论问题不做详细阐述，不清楚的请看相关数学书。目的：能应用MATLAB提供的强大函数进行复杂的方程组、微分、积分等运算。注意：MATLAB数值计算的结果为数值型数据，而不是数学上的解析表达式。线性代数：1.LU分解一个矩阵可以分解为一个上三角矩阵和一个下三角矩阵的乘积，称之为LU分解。LU分解是用高斯主元消去法实现的，通常要对主元位置进行交换，主元交换的方法是将被分解矩阵左乘一个由0-1构成的行交换阵。【调用格式】[L,U,P]=lu(X)对矩阵X进行LU分解，并进行主元交换，[L,U]=lu(X)对矩阵X进行LU分解，无主元交换【说明】L为主对角元素为1的下三角矩阵，U为上三角矩阵，P为行交换矩阵。2.行列式和求逆矩阵的行列式和求逆可以通过LU分解的方法求解。【调用格式】d=det(X)求矩阵X的行列式Y=inv(X)求矩阵X的逆矩阵例1对矩阵进行LU分解、求行列式和求逆操作，请查看计算结果。A=[1,2,3;2,2,3;9,7,5];[L1,U1]=lu(A)；%不带主元交换的LU分解，L1通常不是下三角阵[L2,U2,P]=lu(A)%带主元交换的LU分解，L2为下三角阵det(A)%行列式值Y=inv(A)%矩阵求逆3.特征值和特征向量eig函数用于求解矩阵的特征值和特征向量。【调用格式】D=eig(A)计算矩阵A的特征值，D为特征值构成的向量[V,D]=eig(A)计算矩阵A的特征值对角阵D和特征向量矩阵V[V,D]=eig(A,'nobalance')当矩阵A中有与截断误差近似的数值，用本指令例2计算B的特征值及特征向量，观察结果。B=[-2,1,1;0,2,0;-4,1,3];D=eig(B);[V1,D1]=eig(B,'nobalance');%计算误差为1.0e-014数量级4.线性方程组的解在线性方程组中，独立方程的个数等于独立未知参数的个数，称为恰定方程；独立方程的个数大于独立未知参数的个数，称为超定方程；独立方程的个数小于独立未知参数的个数，称为欠定方程。形如Ax=b的线性方程组可以用以下方法求解：a.左除运算符法左除运算符法的形式为:x=A\b对于一般的非奇异矩阵A，可以求得唯一数值解。欠定方程和超定方程，可以获得最小二乘解。b.广义逆法如果用左除运算符求解的时候出现提示矩阵A为非奇异的警告或者解中出现Nan，则可以采用广义逆法，形式为:x=pinv(A)*bc.符号计算法可以求得方程组的符号解，对于欠定方程可以求得具有自由变量的解。例3求以下3个方程组的解,观察解是否正确。I：10325943213364zyxzyxzyxII：943213364zyxzyxIII：10259321364yxyxyx解：方程I为恰定方程，用左除运算符可以求解A1=[4,6,3;2,3,4;5,2,3];b1=[13;9;10];x1=A1\b1方程II为欠定方程，用3种方法分别求解A2=A1;b2=b1;A2(3,:)=[];b2(3)=[];x2=A2\b2%左除运算符法，求得一组特解，解中非零元素最少x2=01.66671.0000x2=pinv(A2)*b2%广义逆法，求得一组特解x2=0.76921.15381.0000x2=solve('4*x+6*y+3*z=13,2*x+3*y+4*z=9','x,y,z');%符号运算法求通解x2=[x2.x;x2.y;x2.z]x2=5/2-3/2*yy1方程III为超定方程，用左除运算符法求最小二乘解A3=A1;b3=b1;A3(:,3)=[];x3=A3\b3x3=1.45451.3636函数分析1.多项式的根通过roots函数来求取多项式全部的根。【调用格式】r=roots(p)多项式求根函数【说明】p为多项式的系数行向量，r为多项式所有根构成的列向量。2.一元函数零点fzero函数求取一元函数的精确零点。【调用格式】[x,fval,exitflag,output]=fzero(fun,x0,options)一元函数零点的完整调用格式x=fzero(fun,x0)一元函数零点的最简调用格式【说明】1.fzero只能求得x0附近的单个零点，不能求取函数的所有零点。2.输入变量fun表示一元函数，可以是字符串、内联函数或者函数句柄。3.输入变量x0为零点的初始猜测值（自变量值）。如果x0为标量，则求距离x0最近的那个零点；如果x0=[a,b]，要求fun(a)和fun(b)异号，此时求自变量在[a,b]区间内的零点。4.输入变量options为优化迭代选项，是一个结构体。5.输出变量x为零点处的自变量值，输出变量fval为零点处的函数值。6.输出变量exitflag表示函数中止计算的条件。若exitflag0表示找到零点后退出。输出变量output表示程序运行信息，是一个结构体。例1求函数)4.32sin(5.121)(tetft在t0区间内的所有零点。解：运行下列程序，请观察执行结果。y=inline('1-12.5*exp(-t)*sin(2*t+3.4)','t');%构造内联函数t=0:0.1:10;%定义自变量取值区间yv=feval(vectorize(y),t);%计算自变量区间的函数值plot(t,yv);axistight;gridon;[tt,yy]=ginput%请在在函数曲线上用鼠标拾取所有零点，按回车键结束tt=yy=0.00850.0085[t1,fv1,flag]=fzero(y,tt(1));%求第一个拾取点的精确零点[t2,fv2,flag]=fzero(y,tt(2));%求第二个拾取点的精确零点附加学习及练习：内联函数用户可以用M文件来建立函数，函数的功能可以很复杂，函数的输出变量也可以有多个。对于简单的数学表达式，用M文件来建立函数就显得不够方便。MATLAB提供了内联函数的功能，内联函数可以将表达式转换为函数。内联函数是MATLAB面向对象的一个类，其类型名为inline。1.内联函数的建立(1)g=inline('expr')将串表达式expr转换为内联函数(2)g=inline('expr','arg1','arg2',...)将串表达式expr转化为以arg1、arg2等为自变量（输入变量）的内联函数(3)g=inline('expr',n)将串表达式expr转化为以自变量x，P1，P2，…，Pn为自变量的内联函数。其中P必须大写。例阅读下列程序代码，了解inline的使用方法。g=inline('sin(alpha*x)','x','alpha')g=Inlinefunction:g(x,alpha)=sin(alpha*x)g1=inline('a*sin(x(1))*cos(x(2))','a','x')%自变量为向量，函数值为标量g1=Inlinefunction:g1(a,x)=a*sin(x(1))*cos(x(2))g1(1,[pi/4,pi/4])ans=0.5000%自变量为向量，函数值为向量g3=inline('[sin(x(1));cos(x(2));sqrt(x(1).^2+x(2).^2)]')g3=Inlinefunction:g3(x)=[sin(x(1));cos(x(2));sqrt(x(1).^2+x(2).^2)]g3([3,4])ans=0.1411-0.65365.0000g4=inline('sqrt(x^2+P1^2+P2^2)',2)%inline的第三种格式的使用g4=Inlinefunction:g4(x,P1,P2)=sqrt(x^2+P1^2+P2^2)g4(1,3,5)ans=5.9161feval(g4,1,3,5)%内联函数可以通过feval计算求值ans=5.91612.和内联函数有关的函数class(fun)%获取内联函数的数据类型char(fun)%获取内联函数的计算公式字符串argnames(fun)%获取内联函数的输入变量名字vectorize(fun)%使内联函数具有数组运算规则字符串的求值MATLAB提供了字符串求值的函数，利用这些函数，可以用字符串构造MATLAB的函数和命令，并运行这些字符串命令。1.字符串表达式计算【调用格式】y=eval('expression')计算字符串表达式expression[a1,a2,...]=eval('function(b1,b2,...)')计算函数调用的字符串表达式【说明】eval的输入变量只能是字符串例表达式字符串的计算。解：运行下列语句，观察执行结果。clearx=0.4;cmd=['y=x*eye(3),z=sin(x/2)+cos(x)'];eval(cmd);who2.字符串函数计算【调用格式】[y1,y2,...]=feval('function',x1,...,xn)【说明】(1)'function'只能是函数名，不能是表达式字符串。(2)x1、x2等是调用函数'function'的输入变量，即函数的自变量值。(3)y1、y2等是函数的输出变量，即函数的返回值。例feval的使用方法。解：执行下列语句，观察程序的运行结果x=0:2*pi/50:2*pi;feval('plot',x,sin(x),'r',x,cos(x),'b:');%相当于运行plot(x,sin(x),'r',x,cos(x),'b:')附加学习结束3.多元函数的零点fsolve函数求取多元函数的精确零点，但是必须提供零点的大致位置才能进行数值搜索。【调用格式】[x,fval,exitflag,output]=fsolve(fun,x0,options)多元非线性方程求解x=fsolve(fun,x0)多元非线性方程求解的最简格式【说明】1.输入变量fun表示自变量为向量的函数，可以是字符串、内联函数或者函数句柄。2.输入变量x0为零点的初始猜测向量（自变量值）。3.输出变量x为零点处的自变量值，输出变量fval为零点处的函数值。例求方程组0202212211xxexxexx在x1=-5,x2=-5附近的解。解：（1）采用字符串表达函数fun='[2*x(1)-x(2)-exp(-x(1)),-x(1)+2*x(2)-exp(-x(2))]';%字符串方式表达函数x0=[-5;5];[x,fval]=fsolve(fun,x0)（2）采用内联函数表达函数fun=inline('[2*x(1)-x(2)-exp(-x(1)),-x(1)+2*x(2)-exp(-x(2))]','x');x0=[-5;5];[x,fval]=fsolve(fun,x0)4.函数的极值点MATLAB提供了3个求极值点的函数，其输入输出参数的定义和fsolve函数基本相同。【调用格式】[x,fval,exitflag,output]=fminbnd(fun,x1,x2,options)【说明】求一元函数fun在自变量（x1，x2）区间的最小值【调用格式】[x,fval,exitflag,output]=fminsearch(fun,x0,options)【说明】：用单纯形法求多元函数fun在自变量向量x0附近的极小值点【调用格式】[x,fval,exitflag,output]=fminunc(fun,x0,options)【说明】：用拟牛顿法求多元函数fun在自变量向量x0附近的极小值点例求二元函数2)1()(100),(22Xxyyxf在x=-1.2,y=1附近的极小值。解：运行下列程