2019年实验七matlab求解级数有关计算

实验七matlab求解级数有关计算1.级数的基本概念常数项级数：称用加号将数列"“的项连成的式子a+a+ah a+a+ah fah—123为(常数项)无穷级数，简记为n=1。称级数n=1前n项构成的和S=a+a+a+•…+an1二工akk=1收敛，其和为S。limS=S收敛，其和为S。为级数的部分和。若n”n ，则称级数n=1Taylor级数：设函数f(x)在包含x二a的区域内具有各阶导数，则称幕级数£f(n)(a)(x-a)n=f(a)+f'(a)(x-a)+f⑵(a)(x-a)2+•••+f(n)(a)(x-a)n+…n! 2! n!n=0为函数f(X)在x=a的Taylor级数，当a=0时称为Maclaurin(麦克劳林)级数。2.级数的MATLAB命令MATLAB中主要用symsum,taylor求级数的和及进行Taylor展开。symsum(s,v,a,b)表达式s关于变量v从a到b求和taylor(f,a,n)将函数f在a点展为n-1阶Taylor多项式可以用helpsymsum,helptaylor查阅有关这些命令的详细信息例1先用taylor命令观测函数y=sinX的Maclaurin展开式的前几项，例如观测前6项,相应的MATLAB代码为：>>clear;symsx;>>taylor(sin(x),0,1)>>taylor(sin(x),0,2)>>taylor(sin(x),0,3)>>taylor(sin(x),0,4)>>taylor(sin(x),0,5)>>taylor(sin(x),0,6)结果为:ans=0ans=xans=xans=x-1/6*xA3ans=x-1/6*xA3ans=xT/6*x3+l/120*x5然后在同一坐标系里作出函数y二SinX和它的Taylor展开式的前几项构成的多项式函x3 x3x5TOC\o"1-5"\h\zy=x,y=x- ,y=x- + ，•…,数 3! 3! 5! 的图形，观测这些多项式函数的图形向y二sinx的图形的逼近的情况。例如，在区间［0，兀］上作函数y二sinx与多项式函数x3 x3x5y—x,y—x— ,y—x— + -3! 3! 5!图形的MATLAB代码为：>>x=0:0.01:pi;yl=sin(x);y2=x;y3=x-x43/6;y4=x-x43/6+x45/120;>>plot(x,yl,x,y2,':',x,y3,':',x,y4,':')类似地，根据函数的Taylor类似地，根据函数的Taylor级数结果如图3.1，x2 x4 x6cosx—1— + — +•…,xe(—^o,+^o),2! 4! 6!ex=1+x+X2+X3+…,xe(-g,+x),2! 3!x2 x3 x4ln(1+x)—x— + — +•…,xe(—1,1],2 3 4a(a—1)(1+x)a—1+ax+ 2!x2+•…,xe(-1,1).作图观测其展开式的前几项多项式逼近原函数的情况。例2利用幂级数计算指数函数。指数函数可展开为幂级数x2x3 xnex—1+x+ + +•…+—+•…，xe(—g,+8)2! 3! n!其通项为xM/prod(l:n),因此用下列循环相加就可计算出这个级数>>x=input('x=');n=input('n=');y=1;%输入原始数据，初始化y>>fori=1:ny=y+x^i/prod(1:i);end,vpa(y,10),%将通项循环相加，得y执行此程序，分别带入x=1,2,4,-4这四个数，取n=10,y的结果如下2.7l828l80l，7.388994709，54.443l0406，.967l957672e-l而用vpa(exp(1),10),vpa(exp(2),10),vpa(exp(4),10),vpa(exp(-4),10)命令可得"，。‘创‘丘4的⑴位精确有效数字为2.718281828,7.389056099,54.59815003,.1831563889e-1对照可知，用级数法计算的有效数字分别为8，4，2，0位。由此可以看出，这个程序虽然原理上正确，但不好用。对不同的x,精度差别很大。其他存在的问题有：这个程序不能用于x的元素群运算；当x为负数时，它成为交错级数，收敛很慢；此n2程序要做2次乘法，n很大时，乘法次数太多，计算速度很低；对不同的x,要取不同的n才能达到精度要求，因此n不应由用户输入，应该由软件按精度要求来选。正对上面的四个问题，可以采用下面四种方法改进：(1)允许数组输入，改进输出显示x=input('x=');n=input('n=');y=ones(size(x));%输入原始数据，初始化yfori=1:ny=y+x.Ai/prod(1:i);%循环相力口s1=sprintf('%13.0f',i);s2=sprintf('%15.8f',y);%将结果变为字符串disp([s1,s2]) %显示end,执行此程序，输入x=[124-4]，n=10，结果为12.000000003.000000005.00000000-3.0000000022.500000005.0000000013.000000005.0000000032.666666676.3333333323.66666667-5.6666666742.708333337.0000000034.333333335.0000000052.716666677.2666666742.86666667-3.5333333362.718055567.3555555648.555555562.1555555672.718253977.3809523851.80634921-1.0952381082.718278777.3873015953.431746030.5301587392.718281537.3887125254.15414462-0.19223986102.718281807.3889947154.443104060.09671958可以利用exp(-x)=l/exp(x)来避免交错级数的计算；为了减少乘法次数，设一个中间变量z,它的初始值为z=ones(size(x))，把循环体中的计算与句改为y=y+z;z=x.*z/i;这样，求得的z就是z=x.T/i!，于是每个循环只需做一次乘法，计算整个级数只需n次乘法。按这种计算，y的初始值改为y=zeros(size(x))(4)为了按精度选择循环次数，不该使用for循环，而用while语句，它可以设置循环的条件语句，通常可用y+z-y>tol,tol是规定的允许误差.只要相邻的两次y值之差大于tol,循环就继续进行，直到小于tol为止.xexp(x)二(expG-))k当x较大时,exp(x)仍能很快收敛，还可以利用关系式 k，令x1=x/k.k通常取大于x而最接近x的2的幕，例如x=100，就取k=128，可以保证x1的绝对值小于1,这时级数收敛得很快••从练习中可以看出,n取10时(即级数取10项)就能保证7位有效数，而exp(X1)128可以化成X=(…((exp(X1))2)2…)2,即exp(x1)的7次自乘，总共用17次乘法就可完成exp(100)=(…((exp(1°0/128))2)2…)2的计算，这既保证了精度，又提高了速度.例3 编写任意函数展开为各阶泰勒级数的程序，并显示其误差曲线.对于任意函数y=f(x),其泰勒展开式为f(x)=f(a)+f'(a)(x-a)+-凹(x-a)2h f-凹(x-a)n+R(x).2! n! n其中Rn(x)为余项，也就是泰勒展开式的误差.MATLAB语句为>>fxs=input('输入y=f(x)的表达式','s'); %输入原始条件,fxs是字符串>>K=input('输入泰勒级数展开式的阶K');>>a=input('展开的位置a=');>>b=input('展开的区间半宽度b=');>>x=linspace(a-b,a+b);%构成自变量数组,确定其长度和步长>>lx=length(x);dx=2*b/(lx-1);>>y=eval(fxs);%求出y的准确值>>subplot(1,2,1),plot(x,y,'.'),holdon%y的准确值用点线绘出%求出a点的一阶导数，注意求导后数组长度减少1>>Dy=diff(y)/dx;Dya(1)=Dy(round(lx-1)/2);>>yt(1,:)=y(round(lx/2))+Dya(1)*(x-a); %求y的一阶泰勒展开,绘图>>plot(x,yt(1,:))>>fork=2:K>>Dy=diff(y,k)/(dx^k);Dya(k)=Dy(round(lx-k)/2);%求a点k阶导数>>yt(k,:)=yt(k-1,:)+Dya(k)/prod(1:k)*(x-a).^k; %求y的k阶导

>>plot(x,yt(k,:));%绘图>>e(k,:)=y-yt(k,:);%求出yt的误差>>end>>title([fxs,'的各阶泰勒级数曲线']),%注意如何组成标注的字符串>>grid,holdoff,subplot(1,2,2)>>fork=1:Kplot(x,e(k,:)),holdon,end%绘制误差曲线>>title([fxs,'的各阶泰勒级数误差曲线']),grid,holdoff执行此程序，输入fxs=cos(x),K=5,a=0.5,b=2,所得曲线见图3.2(又变为误差曲线).读者可以改变其坐标系范围以仔细观测最关心的部分,也可输入其他函数做验算,注意输入函数应符合元素群运算规则.艺丄.例4 计算级数n=in2的值，可用symsum命令,相应的MATLAB代码为:>>clear;symsk;>>simple(symsum(1/kA2,1,Inf))%simple求解最简形式,Inf为无穷大艺丄n艺丄n4n=1兀490结果为:ans=1/6*piA2类似地可验证艺1兀6n6 945n=1兰1可以猜想有n=1n22,2=1,2,其中m2是正整数，请验证.艺1=e注：可用公式n=0n! 来计算e的近似值。如果要精确到小数点后15位，相应的MATLAB代码为:

>>digits(20); %设置今后数值计算以20位相对精度进行>>a=1.0;kk=1.0; %赋初值>>forn=1:20,kk=kk/n;,a=a+kk;,end>>vpa(a,17) %以17位相对精度给出a的值结果为e沁2.718281828459045.例5(调和级数例5(调和级数)自然数的倒数组成的数列111231n 称为调和数列，由调和数艺1 艺-列构成的级数n=1n称为调和级数，我们把它的前n项部分和k=1k记为H(n)。计算n=12，…,10时H(n)和Inn,ln(n+!)的值，并计算它们的差C(n)=H(n)—lnn，c(n)=H(n)—ln(n+1)，相应的MATLAB代码为:>>H(1)=1;C(1)=1;>>forn=2:100,H(n)=H(n-1)+1/n;,%for为循环语句>>c(n)=H(n)-log(n+1);,C(n)=H(n)-log(n);,end注意观测C(n)单调递减、C(n)单调递增，二者相互接近的现象。计算n=10m(m=1,2,…,6)时C(n)和c(n)的值，注意观测C(n)单调递减、c(n)单调递增，二者趋于同一极限的现象。并求出这个常数C。极限C=lim(1+ +—+•…+——lnn)=0.5771•…ntr 23n称为欧拉(Euler)常数，显然，可以证明它是一个无理数。c(n)=H(n)—ln(n+1)<C(n)=H(n)—lnn,C(n)—C(n)—c(n)=ln(1+丄)n趋于o,故C(n),c(n)趋于同一个极限C习题16-71.用taylor命令观测函数y=f(x)的Maclaurin展开式的前几项，然后在同一坐标系里作出函数y=f(X)和它的Taylor展开式的前几项构成的多项式函数的图形，观测这些多项式函数的图形向y=f(X)的图形的逼近的情况(1)f(x)=arcsinx(2)(1)f(x)=arcsinx(2)f(x)=arctanx(3)f(x)=ex2ex(4)f(x