实验四MATLAB在方程求解和级数中的应用

实验四MATLAB在方程求解和级数中的应用线性映射的迭代与特征向量的计算级数方程和方程组的求解南京邮电大学NanjingUniversityofPostsandTelecommunications一、利用MATLAB进行级数运算的方法和技能在高等数学中，级数一般分为三个部分来叙述，即常数项级数的求和和审敛法则、幂级数的审敛和将函数展开为幂级数、傅立叶级数的性质和将函数展开为傅立叶级数。南京邮电大学NanjingUniversityofPostsandTelecommunications1.常数项级数的求和与审敛在讨论常数项级数时，一般认为，如果级数的部分和的极限存在，则称该级数收敛，并称此极限为级数的和。在MATLAB中，用于级数求和的命令是symsum()，该命令的应用格式为：symsum(comiterm,v,a,b)其中：comiterm为级数的通项表达式，v是通项中的求和变量，a和b分别为求和变量的起点和终点。如果a，b缺省，则v从0变到v-1，如果v也缺省，则系统对comiterm中的默认变量求和。1iia1niia南京邮电大学NanjingUniversityofPostsandTelecommunications例1：求级数，的和。解：利用MATLAB函数symsum设计如下程序：11212nnnI211(21)nInnclearsymsnf1=(2*n-1)/2^n;f2=1/(n*(2*n+1));I1=symsum(f1,n,1,inf)I2=symsum(f2,n,1,inf)运行结果为：I1=3I2=2-2*log(2)南京邮电大学NanjingUniversityofPostsandTelecommunications本例是收敛的情况，如果发散，则求得的和为inf，因此，本方法就可以同时用来解决求和问题和收敛性问题。例2：求级数，的和。解：MATLAB程序如下：141(1)nnnxInclearsymsnxf3=sin(x)/n^2;f4=(-1)^(n-1)*x^n/n;I3=symsum(f3,n,1,inf)I4=symsum(f4,n,1,inf)321sinnxIn南京邮电大学NanjingUniversityofPostsandTelecommunications运行结果为：I3=1/6*sin(x)*pi^2I4=log(1+x)从这个例子可以看出，symsum()这个函数不但可以处理常数项级数，也可以处理函数项级数。南京邮电大学NanjingUniversityofPostsandTelecommunications2.函数的泰勒展开级数是高等数学中函数的一种重要表示形式，有许多复杂的函数都可以用级数简单地来表示，而将一个复杂的函数展开成幂级数并取其前面的若干项来近似表达这个函数是一种很好的近似方法，在学习级数的时候，将一个函数展开成级数有时是比较麻烦的，现在介绍利用MATLAB展开函数的方法。(泰勒级数逼近计算器taylortool)南京邮电大学NanjingUniversityofPostsandTelecommunications在MATLAB中，用于幂级数展开的函数为taylor()，其具体格式为：taylor(function,n,x,a)function是待展开的函数表达式，n为展开项数，缺省时展开至5次幂，即6项，x是function中的变量，a为函数的展开点，缺省为0，即麦克劳林展开。南京邮电大学NanjingUniversityofPostsandTelecommunications例3：将函数sin（x）展开为x的幂级数，分别展开至5次和20次。解：MATLAB程序为：clearsymsxf=sin(x);taylor(f)taylor(f,20)结果为：ans=x-1/6*x^3+1/120*x^5南京邮电大学NanjingUniversityofPostsandTelecommunicationsans=x-1/6*x^3+1/120*x^5-1/5040*x^7+1/362880*x^9-1/39916800*x^11+1/6227020800*x^13-1/1307674368000*x^15+1/355687428096000*x^17-1/121645100408832000*x^19南京邮电大学NanjingUniversityofPostsandTelecommunications例4：将函数(1+x)m展开为x的幂级数，为任意常数。展开至4次幂。解：MATLAB程序为：clearsymsxmf=(1+x)^m;taylor(f,5)运行结果为：ans=1+m*x+1/2*m*(m-1)*x^2+1/6*m*(m-1)*(m-2)*x^3+1/24*m*(m-1)*(m-2)*(m-3)*x^4南京邮电大学NanjingUniversityofPostsandTelecommunications例5：将函数展开为的幂级数。解：MATLAB程序为：clearsymsxf=1/(x^2+5*x-3);taylor(f,5,x,2)pretty(ans)(2)x南京邮电大学NanjingUniversityofPostsandTelecommunications结果为：ans=29/121-9/121*x+70/1331*(x-2)^2-531/14641*(x-2)^3+4009/161051*(x-2)^429702531340094------9/121x+------(x-2)--------(x-2)+----------(x-2)121133114641161051南京邮电大学NanjingUniversityofPostsandTelecommunications南京邮电大学NanjingUniversityofPostsandTelecommunications二、线性方程组、非线性方程、非线性方程组的求解。在MATLAB中，由函数solve()、null()、fsolve()，fzero等来解决线性方程（组）和非线性方程（组）的求解问题，其具体格式如下：X=solve(‘eqn1’,’eqn2’,…,’eqnN’,’var1’,’var2’,…,’varN’)X=fsolve(fun,x0,options)函数solve用来解符号方程、方程组，以及超越方程，如三角函数方程等非线性方程。参数’eqnN’为方程组中的第N个方程，’varN’则是第N个变量。南京邮电大学NanjingUniversityofPostsandTelecommunications函数null(A)则用来解线性方程组AX=O的基础解系，实际是求系数矩阵A的零空间，在null函数中可加入参数’r’，表示有理基。通过求系数矩阵的秩和增广矩阵的秩，可以判定方程组是否有解，以及是否需要求基础解系。另外，还可以用函数fzero来求解非线性方程。用法与fsolve类似。南京邮电大学NanjingUniversityofPostsandTelecommunications例1：求解方程的MATLAB程序为：X=solve(‘x^2-x-6=0’,’x’)结果为：X=3,-2例2：求解方程组的程序为：260xx226060xyyx[X,Y]=solve('x^2+y-6=0','y^2+x-6=0','x','y')结果为：X=2,-3,1/2-1/2*21^(1/2),1/2+1/2*21^(1/2)Y=2,-3,1/2+1/2*21^(1/2),1/2-1/2*21^(1/2)南京邮电大学NanjingUniversityofPostsandTelecommunications例3：求解方程组的程序为：clearformatratA=[5,0,4,2;1,-1,2,1;4,1,2,0;1,1,1,1]；B=[3;1;1;0]；X=A\B134123412312345423214210xxxxxxxxxxxxxx南京邮电大学NanjingUniversityofPostsandTelecommunications例4：求方程组的通解的程序为：clearformatratA=[1,2,2,1;2,1,-2,-2;1,-1,-4,-3]C=null(A,'r')％求出矩阵A的解空间的有理基。结果如下：1234123412342202220430xxxxxxxxxxxx南京邮电大学NanjingUniversityofPostsandTelecommunicationsC=25/3-2-4/31001接着，用命令：symsk1k2X=k1*C(:,1)+k2*C(:,2)南京邮电大学NanjingUniversityofPostsandTelecommunications求出的通解为：X=[2*k1+5/3*k2][-2*k1-4/3*k2][k1][k2]南京邮电大学NanjingUniversityofPostsandTelecommunications例5：求方程组的通解的程序为：clearformatratA=sym('[1,2,2,1;2,1,-2,-2;1,-1,-4,-3]')b=sym('[1;2;2]')B=[A,b]n=length(A(1,:))RA=rank(eval(A))RB=rank(eval(B))1234123412342212222433xxxxxxxxxxxx南京邮电大学NanjingUniversityofPostsandTelecommunicationsif(RA==RB&RA==n)X=eval(A\B)％在方程组满秩时，求出唯一解elseif(RA==RB&RAn)C=eval(A\b)％在方程组不满秩时，求出特解D=null(eval(A),‘r’)％求出矩阵A的零空间的基，即方程组的基础解系symsk1k2X=k1*D(:,1)+k2*D(:,2)+C％求出方程组的全部解elsefprintf('NoSolutionfortheEquations')end南京邮电大学NanjingUniversityofPostsandTelecommunications现在转而来看非线性方程组的求解，对于非线性方程组，用函数fsolve来求解。例6：求解非线性方程组时，采用如下的方法，先建立存放函数的m文件，文件名必须与函数名一致，这里就应该为fun.m，内容如下：1122120.5sin0.3cos00.5cos0.3sin0xxxxxx南京邮电大学NanjingUniversityofPostsandTelecommunicationsfunctiony=fun(x)y(1)=x(1)-0.5*sin(x(1))-0.3*cos(x(2))y(2)=x(2)-0.5*cos(x(1))+0.3*sin(x(2))接着，我们建立另一个m文件fsolve1.m，其内容为：clearformatshortx0=[0.1,0.1]fsolve(@fun,x0,optimset(‘fsolve’))％这里的optimset(‘fsolve’)部分是优化设置，可以不用结果是：0.5414，0.3310。南京邮电大学NanjingUniversityofPostsandTelecommunications三线性映射的迭代与特征向量的计算1.定义关系式(n)mmm(n)2m2(n)1m11)(nm(n)m2m(n)222(n)1211)(n2(n)m1m(n)212(n)1111)(n1xaxaxaxxaxaxaxxaxaxax将向量映射为向量T(n)m(n)2(n)1)x,,x,(xT1