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

1.级数的基本概念

常数项级数：称用加号将数列

an的项连成的式子

a1a2a3ann为（常数项）无穷级数，简记为

an1。称级数

an1n前n项构成的和

nSna1a2a3anak1k

为级数的部分和。若

limSnSn，则称级数

an1n收敛，其和为S。

Taylor级数：设函数f(x)在包含xa的区域内具有各阶导数，则称幂级数

n0f(n)(a)n!(xa)f(a)f'(a)(xa)nf(2)(a)2!(xa)2f(n)(a)n!(xa)n

为函数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多项式 可以用help symsum, help taylor查阅有关这些命令的详细信息 例1 先用taylor命令观测函数ysinx的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 =0 ans =x ans =x

ans =x-1/6*x^3 ans =x-1/6*x^3

1

ans =x-1/6*x^3+1/120*x^5

然后在同一坐标系里作出函数ysinx和它的Taylor展开式的前几项构成的多项式函yx,yxx3数

3!,yxx33!x55!,,的图形，观测这些多项式函数的图形向

ysinx的图形的逼近的情况。例如，在区间[0,]上作函数ysinx与多项式函数

yx,yxx33!,yxx33!x55!图形的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，其中实线表示函数ysinx的图形。

1.510.5001234 图3.1 ysinx的泰勒级数

类似地，根据函数的Taylor级数

6cosx1x22!x44!x6!,x(,),ex1xx22!x33!,x(,),4x)xx2x3ln(123x4,x(1,1],(1x)1x(1)2!x2,x(1,1).

作图观测其展开式的前几项多项式逼近原函数的情况。

例2 利用幂级数计算指数函数。指数函数可展开为幂级数

3nex1xx22!x3!xn!,x(,),

2

其通项为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)命令可得e,e,e,e位精确有效数字为

2.718281828, 7.389056099, 54.59815003, .1831563889e-1

对照可知，用级数法计算的有效数字分别为8，4，2，0位。

由此可以看出，这个程序虽然原理上正确，但不好用。对不同的x,精度差别很大。其他存在的问题有：

这个程序不能用于x的元素群运算；当x为负数时，它成为交错级数，收敛很慢；此

n2244的10

程序要做2次乘法，n很大时，乘法次数太多，计算速度很低；对不同的x,要取不同的n才能达到精度要求，因此n不应由用户输入，应该由软件按精度要求来选。

正对上面的四个问题，可以采用下面四种方法改进： （1）允许数组输入，改进输出显示

x=input('x='); n=input('n='); y=ones(size(x)); %输入原始数据，初始化y

for i=1:n

y=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.00000000 2 2.50000000 5.00000000 13.00000000 5.00000000 3 2.66666667 6.33333333 23.66666667 -5.66666667 4 2.70833333 7.00000000 34.33333333 5.00000000 5 2.71666667 7.26666667 42.86666667 -3.53333333 6 2.71805556 7.35555556 48.55555556 2.15555556 7 2.71825397 7.38095238 51.80634921 -1.09523810 8 2.71827877 7.38730159 53.43174603 0.53015873 9 2.71828153 7.38871252 54.15414462 -0.19223986 10 2.71828180 7.38899471 54.44310406 0.09671958

（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语句,它可以设置循环的条件

3

语句,通常可用y+z-y>tol,tol是规定的允许误差.只要相邻的两次y值之差大于tol,循环就继续进行,直到小于tol为止.

exp(x)(exp(xk))k当x较大时,exp(x)仍能很快收敛,还可以利用关系式,令x1=x/k.k通

常取大于x而最接近x的2的幂,例如x=100,就取k=128,可以保证x1的绝对值小于1,这时级数收敛得很快..从练习中可以看出,n取10时(即级数取10项)就能保证7位有效数,而

exp(x1)128可以化成x(((exp(x1)))),即exp(x1)的7次自乘,总共用17次乘法就可

222222完成exp(100)(((exp(100/128))))的计算,这既保证了精度,又提高了速度. 例3 编写任意函数展开为各阶泰勒级数的程序,并显示其误差曲线.对于任意函数y=f(x),其泰勒展开式为

f(x)f(a)f'(a)(xa)f(2)(a)2!(xa)2f(n)(a)n!(xa)Rn(x).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,'.'), 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(又变为误差曲线).读者可以改变其坐标系范围以仔细观测最关心的部分,也可输入其他函数做验算,注意输入函数应符合元素群运算规则.

4

cos(x)的各阶泰勒级数曲线21.510.5cos(x)的各阶泰勒级数误差曲线1.210.80.600.4-0.5-1-1.5-2-20.20024-0.2-2024 图3.1 ycosx的泰勒级数及误差曲线

例4 计算级数n1>>clear; syms k;

1n2.的值,可用symsum命令,相应的MATLAB代码为：

>>simple(symsum(1/k^2,1,Inf)) %simple求解最简形式,Inf为无穷大

结果为:

ans =1/6*pi^2

类似地可验证

n12k1n4490,n11n66945

.可以猜想有n11n2kmk,k1,2,,其中

emk是正整数,请验证.

注：可用公式n0MATLAB代码为：

1n!来计算e的近似值。如果要精确到小数点后15位，相应的

>>digits(20); %设置今后数值计算以20位相对精度进行 >>a=1.0; kk=1.0; %赋初值

>>for n=1:20, kk=kk/n;, a=a+kk;, end >>vpa(a,17) %以17位相对精度给出a的值

结果为e2.718281828459045.

1,111,,,,23n称为调和数列，由调和数

例5 (调和级数 ) 自然数的倒数组成的数列

5

列构成的级数

n11nn称为调和级数，我们把它的前n项部分和k1k记为H(n)。

1计算n1,2,,100时H(n)和lnn,ln(n!)的值，并计算它们的差

C(n)H(n)lnn，c(n)H(n)ln(n1)，相应的MATLAB代码为：

>>H(1)=1; C(1)=1;

>>for n=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10(m1,2,,6)m时C(n)和c(n)的值，注意观测C(n)单调递减、c(n)单调

递增，二者趋于同一极限的现象。并求出这个常数C。 极限

Clim(1n12131nlnn)0.5771称为欧拉(Euler)常数，可以证明它是一个无理数。

显然，

c(n)H(n)ln(n1)C(n)H(n)lnn, C(n)c(n)ln(11)n趋于0，故C(n),c(n)趋于同一个极限C。

而由于当n时

习题16-7

1. 用taylor命令观测函数yf(x)的Maclaurin展开式的前几项, 然后在同一坐标系里作出函数yf(x)和它的Taylor展开式的前几项构成的多项式函数的图形，观测这些多项式函数的图形向yf(x)的图形的逼近的情况

(1) f(x)arcsinx (2) f(x)arctanx (3) f(x)ef(x)exx2

(4)

f(x)sin2x (5)

1x (6) f(x)ln(x1x)

22. 求公式n11n2k2kmk,k1,2,,中的数

mk,k4,5,6,7,8的值.

3. 利用公式n0n!

1e来计算e的近似值。精确到小数点后100位,这时应计算到这个无

6

穷级数的前多少项?请说明你的理由.

4. 用练习3中所用观测法判断下列级数的敛散性

(1)

n11nnn!n (6)

n23 (2)

n11n2n (3)

n1sin1n (4)

n1lnn1n3

(5)

n1(lnn)n31n (7)

nlnnn11 (8)

n1(1)nn1

2n 7