最新实验七matlab求解级数有关计算

实验七m a t l a b求解级数有关计算实验七 matlab 求解级数有关计算1.级数的基本概念常数项级数：称用加号将数列na 的项连成的式子+++++n a a a a 321为（常数项）无穷级数，简记为∑∞=1n na。

称级数∑∞=1n na前n 项构成的和∑==++++=nk kn n a a a a a S 1321为级数的部分和。

若SS n n =∞→lim ，则称级数∑∞=1n na收敛，其和为S 。

Taylor 级数：设函数)(x f 在包含a x =的区域内具有各阶导数，则称幂级数+-++-+-+=-∑∞=n n n n n a x n a f a x a f a x a f a f a x n a f )(!)()(!2)())((')()(!)()(2)2(0)(为函数)(x f 在a x =的Taylor 级数，当0=a 时称为Maclaurin(麦克劳林)级数。

2．级数的MATLAB 命令MATLAB 中主要用symsum,taylor 求级数的和及进行Taylor 展开。

可以用help symsum, help taylor 查阅有关这些命令的详细信息例1 先用taylor 命令观测函数x y sin =的Maclaurin 展开式的前几项,例如观测前6项, 相应的MATLAB 代码为：>>clear; syms x; >>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 =x ans =xans =x-1/6*x^3 ans =x-1/6*x^3ans =x-1/6*x^3+1/120*x^5然后在同一坐标系里作出函数x y sin =和它的Taylor 展开式的前几项构成的多项式函数,,!5!3,!3,533 x x x y x x y x y +-=-==的图形，观测这些多项式函数的图形向x y sin =的图形的逼近的情况。

例如，在区间],0[π上作函数xy sin =与多项式函数!5!3,!3,533x x x y x x y x y +-=-==图形的MATLAB 代码为： >>x=0:0.01:pi; y1=sin(x); y2=x; y3=x-x.^3/6; y4=x-x.^3/6+ x.^5/120;>>plot(x,y1,x,y2,’:’,x,y3, ’:’,x,y4,’:’)结果如图3.1，其中实线表示函数x y sin =的图形。

图3.1 x y sin =的泰勒级数类似地，根据函数的Taylor 级数).1,1(,!2)1(1)1(],1,1(,432)1ln(),,(,!3!21),,(,!6!4!21cos 243232642-∈+-++=+-∈+-+-=++∞-∞∈++++=+∞-∞∈+-+-=x x x x x x x x x x x x x x e x x x x x xαααα作图观测其展开式的前几项多项式逼近原函数的情况。

例2 利用幂级数计算指数函数。

指数函数可展开为幂级数),,(,!!3!2132+∞-∞∈++++++=x n x x x x e nx其通项为x^n/prod(1:n),因此用下列循环相加就可计算出这个级数>>x=input('x='); n=input('n='); y=1; %输入原始数据，初始化y >>for i=1:n y=y+x^i/prod(1:i); end, vpa(y,10), %将通项循环相加，得y执行此程序，分别带入x=1,2,4,-4这四个数，取n=10,y 的结果如下2.718281801， 7.388994709， 54.44310406， .9671957672e-1而用vpa(exp(1),10), vpa(exp(2),10), vpa(exp(4),10), vpa(exp(-4),10)命令可得442,,,-e e e e 的10位精确有效数字为2.718281828, 7.389056099, 54.59815003, .1831563889e-1对照可知，用级数法计算的有效数字分别为8，4，2，0位。

由此可以看出，这个程序虽然原理上正确，但不好用。

对不同的x,精度差别很大。

其他存在的问题有：这个程序不能用于x 的元素群运算；当x 为负数时，它成为交错级数，收敛很慢；此程序要做22n 次乘法，n 很大时，乘法次数太多，计算速度很低；对不同的x,要取不同的n 才能达到精度要求，因此n 不应由用户输入，应该由软件按精度要求来选。

正对上面的四个问题，可以采用下面四种方法改进： （1）允许数组输入，改进输出显示x=input('x='); n=input('n='); y=ones(size(x)); %输入原始数据，初始化yfor i=1:ny=y+x.^i/prod(1:i); %循环相加s1=sprintf('%13.0f',i); s2=sprintf('%15.8f',y); %将结果变为字符串disp([s1,s2]) %显示end,执行此程序，输入x=[1 2 4 -4]，n=10，结果为1 2.00000000 3.00000000 5.00000000 -3.000000002 2.50000000 5.00000000 13.00000000 5.0000 00003 2.66666667 6.33333333 23.66666667 -5.666666674 2.70833333 7.00000000 34.33333333 5.0000 00005 2.71666667 7.26666667 42.86666667 -3.533333336 2.71805556 7.35555556 48.55555556 2.1555 55567 2.71825397 7.38095238 51.80634921 -1.095238108 2.71827877 7.38730159 53.43174603 0.5301 58739 2.71828153 7.38871252 54.15414462 -0.1922398610 2.71828180 7.38899471 54.44310406 0.0967 1958（2）可以利用exp(-x)=1/exp(x)来避免交错级数的计算；（3）为了减少乘法次数，设一个中间变量z，它的初始值为z=ones(size(x)),把循环体中的计算与句改为y=y+z; z=x.*z/i;这样，求得的z就是z=x.^i/i!,于是每个循环只需做一次乘法，计算整个级数只需n次乘法。

按这种计算，y的初始值改为y=zeros(size(x))(4) 为了按精度选择循环次数,不该使用for循环,而用while语句,它可以设置循环的条件语句,通常可用y+z-y>tol,tol是规定的允许误差.只要相邻的两次y值之差大于tol,循环就继续进行,直到小于tol为止.当x较大时,exp(x)仍能很快收敛,还可以利用关系式kkxx))(exp()exp(,令x1=x/k.k通常取大于x而最接近x的2的幂,例如x=100,就取k=128,可以保证x1的绝对值小于1,这时级数收敛得很快..从练习中可以看出,n 取10时(即级数取10项)就能保证7位有效数,而128)1ex p(x 可以化成222))))1((ex p(( x x =,即exp(x1)的7次自乘,总共用17次乘法就可完成222))))128/100((ex p(()100ex p( =的计算,这既保证了精度,又提高了速度.例3 编写任意函数展开为各阶泰勒级数的程序,并显示其误差曲线.对于任意函数y=f(x),其泰勒展开式为).()(!)()(!2)())((')()()(2)2(x R a x n a f a x a f a x a f a f x f n n n +-++-+-+= 其中)(x R n 为余项,也就是泰勒展开式的误差.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,'.'), hold on %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,:)) >>for k=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, hold off, subplot(1,2,2)>>for k=1:K plot(x,e(k,:)), hold on, end %绘制误差曲线 >>title([fxs,'的各阶泰勒级数误差曲线']),grid,hold off执行此程序,输入fxs=cos(x),K=5,a=0.5,b=2,所得曲线见图3.2(又变为误差曲线).读者可以改变其坐标系范围以仔细观测最关心的部分,也可输入其他函数做验算,注意输入函数应符合元素群运算规则.图3.1 x y cos =的泰勒级数及误差曲线例4 计算级数.112∑∞=n n 的值,可用symsum 命令,相应的MATLAB 代码为：>>clear; syms k;>>simple(symsum(1/k^2,1,Inf)) %simple 求解最简形式,Inf 为无穷大结果为: ans =1/6*pi^2类似地可验证.9451,901616414ππ==∑∑∞=∞=n n n n可以猜想有,,2,1,1212 ==∑∞=k m n kkn kπ其中k m 是正整数,请验证.注：可用公式e n n =∑∞=0!1来计算e 的近似值。