《实验四MATLAB在方程求解和级数中的应用ppt课件.ppt》由会员分享，可在线阅读，更多相关《实验四MATLAB在方程求解和级数中的应用ppt课件.ppt（53页珍藏版）》请在三一办公上搜索。
1、实验四 MATLAB在方程求解和级数中的应用,线性映射的迭代与特征向量的计算,级数,方程和方程组的求解,Nanjing University of Posts and Telecommunications,一、利用MATLAB进行级数运算的方法和技能,在高等数学中，级数一般分为三个部分来叙述，即常数项级数的求和和审敛法则、幂级数的审敛和将函数展开为幂级数、傅立叶级数的性质和将函数展开为傅立叶级数。,Nanjing University of Posts and Telecommunications,1.常数项级数的求和与审敛,在讨论常数项级数时，一般认为，如果级数的部分和 的极限存在，则称该级
2、数收敛，并称此极限为级数的和。在MATLAB中，用于级数求和的命令是symsum()，该命令的应用格式为： symsum(comiterm,v,a,b)其中：comiterm为级数的通项表达式，v是通项中的求和变量，a和b分别为求和变量的起点和终点。如果a，b缺省，则v从0变到v-1，如果v也缺省，则系统对comiterm中的默认变量求和。,Nanjing University of Posts and Telecommunications,例1：求级数 ， 的和。解：利用MATLAB函数symsum设计如下程序：,clearsyms nf1=(2*n-1)/2n;f2=1/(n*(2*n+1
3、);I1=symsum(f1,n,1,inf)I2=symsum(f2,n,1,inf)运行结果为：I1 =3I2 =2-2*log(2),Nanjing University of Posts and Telecommunications,本例是收敛的情况，如果发散，则求得的和为inf，因此，本方法就可以同时用来解决求和问题和收敛性问题。例2：求级数 ， 的和。解：MATLAB程序如下：,clearsyms n xf3=sin(x)/n2;f4=(-1)(n-1)*xn/n;I3=symsum(f3,n,1,inf)I4=symsum(f4,n,1,inf),Nanjing Universi
4、ty of Posts and Telecommunications,运行结果为：I3 =1/6*sin(x)*pi2I4 =log(1+x)从这个例子可以看出，symsum()这个函数不但可以处理常数项级数，也可以处理函数项级数。,Nanjing University of Posts and Telecommunications,2.函数的泰勒展开,级数是高等数学中函数的一种重要表示形式，有许多复杂的函数都可以用级数简单地来表示，而将一个复杂的函数展开成幂级数并取其前面的若干项来近似表达这个函数是一种很好的近似方法，在学习级数的时候，将一个函数展开成级数有时是比较麻烦的，现在介绍利用MAT
5、LAB展开函数的方法。(泰勒级数逼近计算器taylortool),Nanjing University of Posts and Telecommunications,在MATLAB中，用于幂级数展开的函数为taylor()，其具体格式为： taylor(function,n,x,a)function是待展开的函数表达式，n为展开项数，缺省时展开至5次幂，即6项，x是function中的变量，a为函数的展开点，缺省为0，即麦克劳林展开。,Nanjing University of Posts and Telecommunications,例3：将函数sin（x）展开为 x 的幂级数，分别展开至
6、5次和20次。解：MATLAB程序为：clearsyms xf=sin(x);taylor(f)taylor(f,20),结果为：ans =x-1/6*x3+1/120*x5,Nanjing University of Posts and Telecommunications,ans=x-1/6*x3+1/120*x5-1/5040*x7+1/362880*x9-1/39916800*x11+1/6227020800*x13-1/1307674368000*x15+1/355687428096000*x17-1/121645100408832000*x19,Nanjing University
7、 of Posts and Telecommunications,例4：将函数(1+x)m展开为x 的幂级数，为任意常数。展开至4次幂。 解：MATLAB程序为：clearsyms x mf=(1+x)m;taylor(f,5),运行结果为： ans=1+m*x+1/2*m*(m-1)*x2+1/6*m*(m-1)*(m-2)*x3 +1/24*m*(m-1)*(m-2)*(m-3)*x4,Nanjing University of Posts and Telecommunications,例5：将函数 展开为 的幂级数。解：MATLAB程序为：clearsyms xf=1/(x2+5*x-3
8、);taylor(f,5,x,2)pretty(ans),Nanjing University of Posts and Telecommunications,结果为：ans =29/121-9/121*x+70/1331*(x-2)2-531/14641*(x-2)3+4009/161051*(x-2)4 29 70 2 531 3 4009 4 - - 9/121 x + - (x - 2) - - (x - 2) + - (x - 2) 121 1331 14641 161051,Nanjing University of Posts and Telecommunications,Nan
9、jing University of Posts and Telecommunications,二、线性方程组、非线性方程、非线性方程组的求解。,在MATLAB中，由函数solve()、null()、fsolve()，fzero等来解决线性方程（组）和非线性方程（组）的求解问题，其具体格式如下：X=solve(eqn1,eqn2,eqnN,var1,var2,varN)X=fsolve(fun,x0,options)函数solve用来解符号方程、方程组，以及超越方程，如三角函数方程等非线性方程。参数eqnN为方程组中的第N个方程，varN则是第N个变量。,Nanjing University
10、of Posts and Telecommunications,函数null(A)则用来解线性方程组AX=O的基础解系，实际是求系数矩阵A的零空间，在null函数中可加入参数r，表示有理基。通过求系数矩阵的秩和增广矩阵的秩，可以判定方程组是否有解，以及是否需要求基础解系。另外，还可以用函数fzero来求解非线性方程。用法与fsolve类似。,Nanjing University of Posts and Telecommunications,例1：求解方程 的MATLAB程序为：X=solve(x2-x-6=0,x)结果为：X=3, -2例2：求解方程组 的程序为：,X,Y=solve(x2+
11、y-6=0,y2+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),Nanjing University of Posts and Telecommunications,例3：求解方程组 的程序为：clearformat ratA=5, 0, 4, 2;1, -1, 2, 1;4, 1, 2,0;1,1,1,1；B=3;1;1;0；X=AB,Nanjing University of Posts and Telecommunications,例4：
12、求方程组 的通解的程序为：clearformat ratA=1,2,2,1;2,1,-2,-2;1,-1,-4,-3C=null(A,r) 求出矩阵A的解空间的有理基。结果如下：,Nanjing University of Posts and Telecommunications,C = 2 5/3 -2 -4/3 1 0 0 1接着，用命令：syms k1 k2X=k1*C(:,1)+k2*C(:,2),Nanjing University of Posts and Telecommunications,求出的通解为：X = 2*k1+5/3*k2 -2*k1-4/3*k2 k1 k2,Na
13、njing University of Posts and Telecommunications,例5：求方程组 的通解的程序为：clearformat ratA=sym(1,2,2,1;2,1,-2,-2;1,-1,-4,-3)b=sym(1;2;2)B=A,bn=length(A(1,:)RA=rank(eval(A)RB=rank(eval(B),Nanjing University of Posts and Telecommunications,if(RA=RB&RA=n) X=eval(AB) 在方程组满秩时，求出唯一解elseif(RA=RB&RAn) C=eval(Ab) 在方程
14、组不满秩时，求出特解D=null(eval(A),r) 求出矩阵A的零空间的 基，即方程组的基础解系 syms k1 k2 X=k1*D(:,1)+k2*D(:,2)+C 求出方程组的全部解else fprintf(No Solution for the Equations)end,Nanjing University of Posts and Telecommunications,现在转而来看非线性方程组的求解，对于非线性方程组，用函数fsolve来求解。 例6：求解非线性方程组 时，采用如下的方法，先建立存放函数的m文件，文件名必须与函数名一致，这里就应该为fun.m，内容如下：,Nanj
15、ing University of Posts and Telecommunications,function y=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，其内容为：clearformat shortx0=0.1,0.1fsolve(fun,x0,optimset(fsolve)这里的optimset(fsolve)部分是优化设置，可以不用结果是：0.5414，0.3310。,Nanjing University of Posts a
16、nd Telecommunications,三 线性映射的迭代与特征向量的计算,Nanjing University of Posts and Telecommunications,写成矩阵形式,其中 分别为 与,，A为mm矩阵,形如 y=Ax 的映射称为线性映射.给出一个初始向量 ,将上述映射反复作用可得序列： ， ， ， ， 我们将这一过程称为线性映射的迭代，其中矩阵A称为迭代矩阵 。,Nanjing University of Posts and Telecommunications,2.天气问题,问题1 某地区的天气可分为两种状态：晴、阴雨.若今天的天气为晴，则明天晴的概率为3/4，阴
17、雨的概率为1/4；如果今天为阴雨天，则明天晴的概率为7/18，阴雨的概率为11/18.我们可以用一矩阵来表示这种变化，矩阵称为转移矩阵（这些概率可以通过观察该地区以往几年每天天气变化的测量数据来确定）,试根据这些数据来判断该地区的天气变化情况,Nanjing University of Posts and Telecommunications,设某天是晴的概率为 ，阴雨的概率为 ，这一天的天气状态用向量 来表示，k天之后的天气状态用向量 来表示.,则由全概率公式可以得到：,即,可得,Nanjing University of Posts and Telecommunications,下面我们来
18、求解，设,A1=3/4,7/18;1/4,11/18;p=0.5;0.5;for i=1:20 p(:,i+1)=A1*p(:,i);endp,得到p =,Nanjing University of Posts and Telecommunications,Columns 1 through 6 0.5000 0.5694 0.5945 0.6036 0.6068 0.6080 0.5000 0.4306 0.4055 0.3964 0.3932 0.3920 Columns 7 through 12 0.6085 0.6086 0.6087 0.6087 0.6087 0.6087 0.39
19、15 0.3914 0.3913 0.3913 0.3913 0.3913 Columns 13 through 18 0.6087 0.6087 0.6087 0.6087 0.6087 0.6087 0.3913 0.3913 0.3913 0.3913 0.3913 0.3913 Columns 19 through 21 0.6087 0.6087 0.6087 0.3913 0.3913 0.3913,Nanjing University of Posts and Telecommunications,第8天之后，晴、阴雨的概率便稳定下来， 均约等于,问题：这个稳定值是否与p0有关？
20、是否与A1有关？,由于 ，因此有 .,由线性代数的知识我们知道，如存在等式 说明A1有一特征值1，而 恰是其对应的特征向量。,Nanjing University of Posts and Telecommunications,命令P,D=eig(A1) 可求得p = 0.8412 -0.7071 0.5408 0.7071,d = 1.0000 0 0 0.3611,说明A1确有一特征值是1，而对应的特征向量是(0.8412,0.5408)T ,与算得的 差距较大。,思考：如何解释这种差异？,实际上，通过计算可得 (0.6087,0.3913)T=0.7236(0.8412,0.5408)T
21、,Nanjing University of Posts and Telecommunications,能不能直接计算出 p(k) 呢？,如果A1有两个不同的特征值（前面计算已知），这个问题容易解决（见板书）,由命令 P,D=eig(sym(A1)可求得A1的特征向量与特征值的精确值，解方程组 可得u , v，从而得到 p (k) 。,Nanjing University of Posts and Telecommunications,3.平面线性映射迭代 (m=2的情况),取 ，,我们编程观察点列 构成的图形,A=4,2;1,3; x=1;2;t=;for i=1:20 x=A*x; t(i
22、,1:2)=x;endplot(t(1:20,1),t(1:20,2),*),Nanjing University of Posts and Telecommunications,由图形可以看出，迭代序列不收敛，但点列似乎呈线性分布。,现在求序列的一般项,由命令P,D=eig(sym(A)，得到,P = -1, 2 1, 1,D = 2, 0 0, 5,于是AP=PD，或A=PDP-1,Nanjing University of Posts and Telecommunications,最后的计算留着作业,Nanjing University of Posts and Telecommunic
23、ations,4.线性迭代的极限性质,对于前例，随机任取一些初始向量，将迭代得到的所有向量画在一张图中,A=4,2;1,3;t=;for i=1:20 x=2*rand(2,1)-1; t(length(t)+1,1:2)=x; for j=1:40 x=A*x; t(length(t)+1,1:2)=x; endendplot(t(:,1),t(:,2),*)grid(on),Nanjing University of Posts and Telecommunications,这些点在一条直线上？如果是这样，那么它们的两个分量的比值应该相等或差异不大，请验证！,Nanjing Univers
24、ity of Posts and Telecommunications,4.线性迭代的归一化,线性迭代序列不一定收敛，由前面讨论知道，迭代序列不收敛时，序列分量似乎有趋于无穷大的倾向。如果我们关心的仅仅是序列每个分量绝对值的大小关系，而不是它们的具体数值，此时每个分量可以同时除以某一个数。,如果每次除以绝对值最大的那个分量，以保证绝对值最大的分量等于1，这个过程称为归一化。,这个过程便是：,Nanjing University of Posts and Telecommunications,.,.,Nanjing University of Posts and Telecommunicatio
25、ns,这样得到的序列 yn 是否收敛呢？,如果 yn 的极限存在，那么m(xn)的极限也存在，(知道为什么吗？)设它们极限分别为y与，容易推导出,即,说明是A的特征值，而y 是对应的特征向量。,Nanjing University of Posts and Telecommunications,定理9.1 设m阶实方阵A有m个线性无关的特征向量 ，A的m个特征值满足下列关系：则对任意的非零初始向量 ,按上述迭代过程得到 及 ，有： （其中a是一个非零常数），,Nanjing University of Posts and Telecommunications,比赛名次问题,六名选手A、B、C、
26、D、E和F进行围棋单循环比赛，其比赛结果如下：A战胜B、D、E、F；B战胜D、E、F；C战胜A、B、D；D战胜E、F；E战胜C、F；F战胜C.请你给六名选手排一个合理的名次.,Nanjing University of Posts and Telecommunications,按通常的方法可以给选手们排一个初步的名次.即每战胜一名选手得1分，按总得分排名次.六名选手的得分排成一个向量为：,我们得到的初步名次为：1：A；2、3：B、C；4、5：D、E；6：F.,那么，这样排名次是否合理？并列名次怎么处理？,Nanjing University of Posts and Telecommunica
27、tions,我们排名次的一个重要根据是每战胜一名选手得1分.事实上，这一根据并不十分合理，胜“强者”的得分比胜“弱者”的得分应该高.那么，谁是“强者”， 谁是“弱者”呢？得分又应该给多少呢？,Nanjing University of Posts and Telecommunications,我们可以按初步得到的得分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),将上述过程进行
28、下去，编程计算sk, k=3,4,5,如何给出各选手的最终名次？,一个直观的想法是计算 sk 的极限，如果极限存在，我们可以依照极限来排名次.,Nanjing University of Posts and Telecommunications,按照前面的方法来计算sk很麻烦，我们想办法来给出 sk 的通式.我们把各选手的成绩写成下面的矩阵形式：,显然Sk=MSk-1 , s0=(1,1,1,1,1,1)T,Nanjing University of Posts and Telecommunications,编程计算Sk,观察Sk的极限是否存在。,由于我们关心的是各选手的名次，向量Sk的各个分
29、量同时除以一个正数并不影响他们的名次排列.因此可以用2.2节中的归一化迭代的方法来求出选手的名次.,由于归一化迭代方法中，序列收敛于迭代矩阵的最大特征值的特征向量，因此在一定条件下我们可以用该特征向量来给出比赛选手的名次,Nanjing University of Posts and Telecommunications,求M的最大特征值的对应的特征向量，用此向量来确定选手的名次.,Nanjing University of Posts and Telecommunications,format short; A=0,1,0,1,1,1; %输入待求的矩阵A 0,0,0,1,1,1;1,1,0
30、,1,0,0;0,0,0,0,1,1;0,0,1,0,0,1;0,0,1,0,0,0; v,d=eigs(A); tbmax=max(d(:); %最大特征值 m,n=size(v); %得到行数和列数,Nanjing University of Posts and Telecommunications,sum = 0; %将特征向量标准化 for i=1:m sum = sum + v(i,1); end tbvector = v(:,1); for i=1:m tbvector(i,1)= v(i,1)/sum; end disp(=); disp(输入的矩阵为：);,Nanjing University of Posts and Telecommunications,A disp(所有的特征向量和特征值为：); v d disp(最大的特征值为：); tbmax disp(最大的特征值对应的特征向量为（标准化后的）：); tbvector,
