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

实验四MATLAB在方程求解

和级数中的应用

线性映射的迭代与特征向量的计算级数方程和方程组的求解实验四MATLAB在方程求解

和级数中的应用线性映射的一、利用MATLAB进行级数运算的方法和技能在高等数学中，级数一般分为三个部分来叙述，即常数项级数的求和和审敛法则、幂级数的审敛和将函数展开为幂级数、傅立叶级数的性质和将函数展开为傅立叶级数。NanjingUniversityofPostsandTelecommunications一、利用MATLAB进行级数运算的方法和技能在高等数学中，级21.常数项级数的求和与审敛在讨论常数项级数时，一般认为，如果级数的部分和的极限存在，则称该级数收敛，并称此极限为级数的和。在MATLAB中，用于级数求和的命令是symsum()，该命令的应用格式为：

symsum(comiterm,v,a,b)其中：comiterm为级数的通项表达式，v是通项中的求和变量，a和b分别为求和变量的起点和终点。如果a，b缺省，则v从0变到v-1，如果v也缺省，则系统对comiterm中的默认变量求和。NanjingUniversityofPostsandTelecommunications1.常数项级数的求和与审敛在讨论常数项级数时，一般认为，如果3例1：求级数，的和。解：利用MATLAB函数symsum设计如下程序：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例1：求级数，4本例是收敛的情况，如果发散，则求得的和为inf，因此，本方法就可以同时用来解决求和问题和收敛性问题。例2：求级数，的和。解：MATLAB程序如下：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本例是收敛的情况，如果发散，则求得的和为inf，因此，本方法5运行结果为：I3=1/6*sin(x)*pi^2I4=log(1+x)从这个例子可以看出，symsum()这个函数不但可以处理常数项级数，也可以处理函数项级数。NanjingUniversityofPostsandTelecommunications运行结果为：NanjingUniversityofPo62.函数的泰勒展开级数是高等数学中函数的一种重要表示形式，有许多复杂的函数都可以用级数简单地来表示，而将一个复杂的函数展开成幂级数并取其前面的若干项来近似表达这个函数是一种很好的近似方法，在学习级数的时候，将一个函数展开成级数有时是比较麻烦的，现在介绍利用MATLAB展开函数的方法。(泰勒级数逼近计算器taylortool)

NanjingUniversityofPostsandTelecommunications2.函数的泰勒展开级数是高等数学中函数的一种重要表7在MATLAB中，用于幂级数展开的函数为taylor()，其具体格式为：

taylor(function,n,x,a)function是待展开的函数表达式，n为展开项数，缺省时展开至5次幂，即6项，x是function中的变量，a为函数的展开点，缺省为0，即麦克劳林展开。NanjingUniversityofPostsandTelecommunications在MATLAB中，用于幂级数展开的函数为taylor()，其8例3：将函数sin（x）展开为x的幂级数，分别展开至5次和20次。解：MATLAB程序为：clearsymsxf=sin(x);taylor(f)taylor(f,20)结果为：ans=x-1/6*x^3+1/120*x^5NanjingUniversityofPostsandTelecommunications例3：将函数sin（x）展开为x的幂级数，分别展开至5次9ans=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^19NanjingUniversityofPostsandTelecommunicationsans=x-1/6*x^3+1/120*x^5-Nanjin10例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^4NanjingUniversityofPostsandTelecommunications例4：将函数(1+x)m展开为x的幂级数，为任意常数。展开11例5：将函数展开为的幂级数。解：MATLAB程序为：clearsymsxf=1/(x^2+5*x-3);taylor(f,5,x,2)pretty(ans)NanjingUniversityofPostsandTelecommunications例5：将函数展开为12结果为：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)121133114641161051NanjingUniversityofPostsandTelecommunications结果为：NanjingUniversityofPost13NanjingUniversityofPostsandTelecommunicationsNanjingUniversityofPostsan14二、线性方程组、非线性方程、非线性方程组的求解。在MATLAB中，由函数solve()、null()、fsolve()，fzero等来解决线性方程（组）和非线性方程（组）的求解问题，其具体格式如下：X=solve(‘eqn1’,’eqn2’,…,’eqnN’,’var1’,’var2’,…,’varN’)X=fsolve(fun,x0,options)函数solve用来解符号方程、方程组，以及超越方程，如三角函数方程等非线性方程。参数’eqnN’为方程组中的第N个方程，’varN’则是第N个变量。NanjingUniversityofPostsandTelecommunications二、线性方程组、非线性方程、非线性方程组的求解。在MATLA15函数null(A)则用来解线性方程组AX=O的基础解系，实际是求系数矩阵A的零空间，在null函数中可加入参数’r’，表示有理基。通过求系数矩阵的秩和增广矩阵的秩，可以判定方程组是否有解，以及是否需要求基础解系。另外，还可以用函数fzero来求解非线性方程。用法与fsolve类似。

NanjingUniversityofPostsandTelecommunications函数null(A)则用来解线性方程组AX=O的基础解系，实际16例1：求解方程的MATLAB程序为：X=solve(‘x^2-x-6=0’,’x’)结果为：X=3,-2例2：求解方程组的程序为：[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例1：求解方程的MATLA17例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\BNanjingUniversityofPostsandTelecommunicationsNanjingUniversityofPostsan18例4：求方程组

的通解的程序为：clearformatratA=[1,2,2,1;2,1,-2,-2;1,-1,-4,-3]C=null(A,'r')％求出矩阵A的解空间的有理基。结果如下：

NanjingUniversityofPostsandTelecommunicationsNanjingUniversityofPostsan19C=25/3-2-4/31001接着，用命令：symsk1k2X=k1*C(:,1)+k2*C(:,2)NanjingUniversityofPostsandTelecommunicationsC=NanjingUniversityofPosts20求出的通解为：X=[2*k1+5/3*k2][-2*k1-4/3*k2][k1][k2]NanjingUniversityofPostsandTelecommunications求出的通解为：NanjingUniversityofP21例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))NanjingUniversityofPostsandTelecommunicationsNanjingUniversityofPostsan22if(RA==RB&RA==n)X=eval(A\B)％在方程组满秩时，求出唯一解elseif(RA==RB&RA<n)C=eval(A\b)％在方程组不满秩时，求出特解D=null(eval(A),‘r’)％求出矩阵A的零空间的基，即方程组的基础解系symsk1k2X=k1*D(:,1)+k2*D(:,2)+C％求出方程组的全部解elsefprintf('NoSolutionfortheEquations')endNanjingUniversityofPostsandTelecommunicationsNanjingUniversityofPostsan23现在转而来看非线性方程组的求解，对于非线性方程组，用函数fsolve来求解。例6：求解非线性方程组时，采用如下的方法，先建立存放函数的m文件，文件名必须与函数名一致，这里就应该为fun.m，内容如下：NanjingUniversityofPostsandTelecommunicationsNanjingUniversityofPostsan24functiony=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。

NanjingUniversityofPostsandTelecommunicationsNanjingUniversityofPostsan25三线性映射的迭代与特征向量的计算1.定义关系式将向量映射为向量NanjingUniversityofPostsandTelecommunications三线性映射的迭代与特征向量的计算1.定义关系26写成矩阵形式其中分别为与，A为m×m矩阵形如y=Ax的映射称为线性映射.给出一个初始向量,将上述映射反复作用可得序列：，，，…，，…我们将这一过程称为线性映射的迭代，其中矩阵A称为迭代矩阵。NanjingUniversityofPostsandTelecommunications写成矩阵形式其中分别为272.天气问题问题1某地区的天气可分为两种状态：晴、阴雨.若今天的天气为晴，则明天晴的概率为3/4，阴雨的概率为1/4；如果今天为阴雨天，则明天晴的概率为7/18，阴雨的概率为11/18.我们可以用一矩阵来表示这种变化，矩阵称为转移矩阵（这些概率可以通过观察该地区以往几年每天天气变化的测量数据来确定）试根据这些数据来判断该地区的天气变化情况NanjingUniversityofPostsandTelecommunications2.天气问题问题1某地区的天气可分为两种状态：晴、阴雨.28设某天是晴的概率为，阴雨的概率为，这一天的天气状态用向量来表示，k天之后的天气状态用向量来表示.则由全概率公式可以得到：即可得NanjingUniversityofPostsandTelecommunications设某天是晴的概率为，阴雨的概率为29下面我们来求解，设

A1=[3/4,7/18;1/4,11/18];p=[0.5;0.5];fori=1:20p(:,i+1)=A1*p(:,i);endp得到p=NanjingUniversityofPostsandTelecommunications下面我们来求解，设A1=[3/4,7/18;1/4,11/30

Columns1through60.50000.56940.59450.60360.60680.60800.50000.43060.40550.39640.39320.3920Columns7through120.60850.60860.60870.60870.60870.60870.39150.39140.39130.39130.39130.3913Columns13through180.60870.60870.60870.60870.60870.60870.39130.39130.39130.39130.39130.3913Columns19through210.60870.60870.60870.39130.39130.3913NanjingUniversityofPostsandTelecommunicationsColumns1through6Nanjing31第8天之后，晴、阴雨的概率便稳定下来，均约等于问题：这个稳定值是否与p0有关？是否与A1有关？由于，因此有.由线性代数的知识我们知道，如存在等式说明A1有一特征值1，而恰是其对应的特征向量。NanjingUniversityofPostsandTelecommunications第8天之后，晴、阴雨的概率便稳定下来，均约等于32命令[P,D]=eig(A1)可求得p=0.8412-0.70710.54080.7071d=1.0000000.3611说明A1确有一特征值是1，而对应的特征向量是(0.8412,0.5408)T,与算得的差距较大。思考：如何解释这种差异？实际上，通过计算可得(0.6087,0.3913)T=0.7236(0.8412,0.5408)TNanjingUniversityofPostsandTelecommunications命令[P,D]=eig(A1)可求得d=说明A1确有一特33能不能直接计算出p(k)呢？如果A1有两个不同的特征值（前面计算已知），这个问题容易解决（见板书）由命令[P,D]=eig(sym(A1))可求得A1的特征向量与特征值的精确值，解方程组可得u,v，从而得到p(k)。NanjingUniversityofPostsandTelecommunications能不能直接计算出p(k)呢？如果A1有两个不同的特征值（343.平面线性映射迭代(m=2的情况)取，我们编程观察点列构成的图形A=[4,2;1,3];x=[1;2];t=[];fori=1:20x=A*x;t(i,1:2)=x;endplot(t(1:20,1),t(1:20,2),’*’)NanjingUniversityofPostsandTelecommunications3.平面线性映射迭代(m=2的情况)取35由图形可以看出，迭代序列不收敛，但点列似乎呈线性分布。现在求序列的一般项由命令[P,D]=eig(sym(A))，得到P=[-1,2][1,1]D=[2,0][0,5]于是AP=PD，或A=PDP-1NanjingUniversityofPostsandTelecommunications由图形可以看出，迭代序列不收敛，但点列似乎呈线性分布。现在求36最后的计算留着作业NanjingUniversityofPostsandTelecommunications最后的计算留着作业NanjingUniversityof374.线性迭代的极限性质对于前例，随机任取一些初始向量，将迭代得到的所有向量画在一张图中A=[4,2;1,3];t=[];fori=1:20x=2*rand(2,1)-1;t(length(t)+1,1:2)=x;forj=1:40x=A*x;t(length(t)+1,1:2)=x;endendplot(t(:,1),t(:,2),'*')grid('on')NanjingUniversityofPostsandTelecommunications4.线性迭代的极限性质对于前例，随机任取一些初始向量，将迭代38这些点在一条直线上？如果是这样，那么它们的两个分量的比值应该相等或差异不大，请验证！NanjingUniversityofPostsandTelecommunications这些点在一条直线上？如果是这样，那么它们的两个分量的比值应该394.线性迭代的归一化线性迭代序列不一定收敛，由前面讨论知道，迭代序列不收敛时，序列分量似乎有趋于无穷大的倾向。如果我们关心的仅仅是序列每个分量绝对值的大小关系，而不是它们的具体数值，此时每个分量可以同时除以某一个数。如果每次除以绝对值最大的那个分量，以保证绝对值最大的分量等于1，这个过程称为归一化。这个过程便是：NanjingUniversityofPostsandTelecommunications4.线性迭代的归一化线性迭代序列不一定收敛，由前面讨论知道，40…………..…………..NanjingUniversityofPostsandTelecommunications…………..…………..NanjingUniversity41这样得到的序列{yn}是否收敛呢？如果{yn

}的极限存在，那么m(xn)的极限也存在，(知道为什么吗？)设它们极限分别为y与，容易推导出即说明是A的特征值，而y是对应的特征向量。NanjingUniversityofPostsandTelecommunications这样得到的序列{yn}是否收敛呢？如果{yn}42定理 9.1设m阶实方阵A有m个线性无关的特征向量，A的m个特征值满足下列关系：则对任意的非零初始向量,按上述迭代过程得到及，有：（其中a是一个非零常数），NanjingUniversityofPostsandTelecommunications定理 9.1设m阶实方阵A有m个线性无关的特征向Nanj43比赛名次问题六名选手A、B、C、D、E和F进行围棋单循环比赛，其比赛结果如下：A战胜B、D、E、F；B战胜D、E、F；C战胜A、B、D；D战胜E、F；E战胜C、F；F战胜C.请你给六名选手排一个合理的名次.NanjingUniversityofPostsandTelecommunications比赛名次问题六名选手A、B、C、D、E和F进行围棋单循环比44按通常的方法可以给选手们排一个初步的名次.即每战胜一名选手得1分，按总得分排名次.六名选手的得分排成一个向量为：我们得到的初步名次为：1：A；2、3：B、C；4、5：D、E；6：F.

那么，这样排名次是否合理？并列名次怎么处理？NanjingUniversityofPostsandTelecommunications按通常的方法可以给选手们排一个初步的名次.即每我们得到的初步45我们排名次的一个重要根据是每战胜一名选手得1分.事实上，这一根据并不十分合理，胜“强者”的得分比胜“弱者”的得分应该高.那么，谁是“强者”，谁是“弱者”呢？得分又应该给多少呢？NanjingUniversityofPostsandTelecommunications我们排名次的一个重要根据是每战胜一名NanjingUniv46我们可以按初步得到的得分S1,判定强弱，得分就按照s1的分量大小给定.比如，A胜了B、D、E、F，而s1中对应于B、D、E、F的分量分别为3,2,2,1.所以A的得分分别为3,2,2,1A的总得分为.类似地可得出其它选手的得分，于是得到得分向量s2=(8,5,9,3,4,3)将上述过程进行下去，编程计算sk,k=3,4,5,…如何给出各选手的最终名次？一个直观的想法是计算sk的极限，如果极限存在，我们可以依照极限来排名次.NanjingUniversityofPostsandTelecommunications我们可以按初步得到的得分S1,判定强弱，得分就按照s1的分47按照前面