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

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

1.级数的基本概念

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

an的项连成的式子

a1?a2?a3???an??

为（常数项）无穷级数，简记为

?an?1?n。称级数

?an?1?n前n项构成的和

nSn?a1?a2?a3???an??akk?1

为级数的部分和。若n??limSn?S，则称级数n?1?a?n收敛，其和为S。

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

?n?0?f(n)(a)f(2)(a)f(n)(a)n2(x?a)?f(a)?f'(a)(x?a)?(x?a)???(x?a)n??n!2!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

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

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

x3x3x5y?x,y?x?,y?x??,?,3!3!5!数的图形，观测这些多项式函数的图形向y?sinx的图形的逼近的情况。例如，在区间[0,?]上作函数y?sinx与多项式函数

x3x3x5y?x,y?x?,y?x??3!3!5!图形的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级数

x2x4x6cosx?1?????,x?(??,??),2!4!6!x2x3xe?1?x????,x?(??,??),2!3!x2x3x4ln(1?x)?x?????,x?(?1,1],234?(??1)2(1?x)??1??x?x??,x?(?1,1).2!

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

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

x2x3xne?1?x???????,x?(??,??),2!3!n!

x其通项为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的10位精确有效数字为

2.718281828, 7.389056099, 54.59815003, .1831563889e-1

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

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

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

24?4n2程序要做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('.0f',i); s2=sprintf('.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.19223

986

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语句,它可以设置循环的条件语句,通常可用y+z-y>tol,tol是规定的允许误差.只要相邻的两次y值之差大于tol,循环就继续进行,直到小于tol为止.

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

常取大于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次乘法就可

222exp(100)?(?((exp(100/128)))?)完成的计算,这既保证了精度,又提高了速度.

例3 编写任意函数展开为各阶泰勒级数的程序,并显示其误差曲线.对于任意函数y=f(x),其泰勒展开式为

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

cos(x)的各阶泰勒级数曲线21.510.50.600.4-0.5-1-1.5-2-20.2cos(x)的各阶泰勒级数误差曲线1.210.80024-0.2-2024 图3.1 y?cosx的泰勒级数及误差曲线

1.?2例4 计算级数n?1n的值,可用symsum命令,相应的MATLAB代码为：

>>clear; syms k;

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

?结果为:

ans =1/6*pi^2

类似地可验证

1?4?,?490nn?1??1?6?.?6945nn?1

?1?2k?,k?1,2,?,?2kmmk可以猜想有n?1n其中k是正整数,请验证.

1?e?n!注：可用公式n?0来计算e的近似值。如果要精确到小数点后15位，相应的

MATLAB代码为：

?

>>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.

1111,,,?,,?n例5 (调和级数 ) 自然数的倒数组成的数列23称为调和数列，由调和数

n11??nnn?1列构成的级数称为调和级数，我们把它的前项部分和k?1k记为H(n)。

?计算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)单调递增，二者相互接近的现象。

mn?10(m?1,2,?,6)时C(n)和c(n)的值，注意观测C(n)单调递减、c(n)单调 计算

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

C?lim(1?n??111?????lnn)?0.5771?23n

称为欧拉(Euler)常数，可以证明它是一个无理数。

显然，

c(n)?H(n)?ln(n?1)?C(n)?H(n)?lnn,

1C(n)?c(n)?ln(1?)n趋于0，故C(n),c(n)趋于同一个极限C。 而由于当n??时

习题16-7

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

xf(x)?arcsinxf(x)?arctanxf(x)?e(1) (2) (3)

2exf(x)?22f(x)?sinxf(x)?ln(x?1?x) 1?x(4) (5) (6)

1?2k?,k?1,2,?,?2km,k?4,5,6,7,8的值. mk2. 求公式n?1n中的数k?1?e?n!3. 利用公式n?0来计算e的近似值。精确到小数点后100位,这时应计算到这个无

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

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

?1?23(1) n?1n?n (2) n!?n(5) n?1n (6)

??1?nn?1n2 (3)

??1sin?n (4) n?1?lnn1?3n?1n

?1?nn?3(lnn) (7)

?1?n?1nlnn (8)

(?1)nn?2n?1n?1

?

本文来源：https://www.bwwdw.com/article/2h15.html