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

1、实验七mat lab求解级数有关计算1. 级数的基本概念常数项级数：称用加号将数列心的项连成的式子6/| +a2+a3+- + an +-8 00为（常数项）无穷级数，简记为称级数前项构成的和S =6 +。2 +3 + + =你A-1wlimS=S工 5c为级数的部分和。若y,则称级数一|收敛，其和为S。Taylor级数：设函数/（）在包含x的区域内具有各阶导数，则称幕级数）”=/）+ 广（）（）+ （4）2 + -+2（1._律 + .缶川2!川为函数/（X）在尤=的Taylor级数，当。=时称为MaclaurinC麦克劳林）级数。2. 级数的MATLAB命令MATLAB中主要用symsum

2、, 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命令观测函数，= 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)ta

3、ylor (sin(x), 0, 5)taylor (sin(x), 0, 6)结果为：ans =0ans =xans =xans =x-l/6*x3ans =x-l/6*x3然后在同一坐标系里作出函数? = sinX和它的Taylor展开式的前几项构成的多项式函X3牙 35y = x,y = x-.y = -v- +数3!3!5! 的图形，观测这些多项式函数的图形向y = sinx的图形的逼近的情况。例如，在区间【，刃上作函数 = sinx与多项式函数X5牙X5y = x.y = x- 9y = x- + 3! 3!5!图形的MATLAB代码为：x=0:0. 01:pi;yl=sin(x)

4、;y2=x; y3二x-x. 3/6; y4二x-x. 3/6+x.5/120;plot (x, yl, x, y2, ： , x, y3, ： , x, y4,：)结果如图3.1,其中实线表示函数 = sinx的图形。图3.1) = sinx的泰勒级数类似地，根据函数的Taylor级数X2 X4X6cosx = 1F+ , x e (一s,+s),2!4!6!x2 疋ex = 1 + x + + + u (s,+s),2!3!34ln(l + x) = x-+ 一 + ，xu (-1,11,234(1 + x)a = 1 + OV+ .v2 +, x e (一1,1) 2!作图观测其展开式的

5、前几项多项式逼近原函数的情况。例2利用幕级数讣算指数函数。指数函数可展开为幕级数X2 X3X”ex = 1 + x + + + + + ，x w (yo,+r),2!3!n!其通项为xn/prod(l:n),因此用下列循环相加就可讣算出这个级数x=input C x=，); n=input C n=,); y=l;玄输入原始数据，初始化 y for i=l:n y=y+x*i/prod(l:i); end, vpa(y, 10), %将通项循环相加,得 y执行此程序，分别带入X二1,2, 4,-4这四个数，取n=10, y的结果如下2.71828180b7.388994709,54.44310

6、406, 9671957672e-l而用 vpa(exp(l), 10), vpa(exp(2), 10), vpa(exp(4), 10), vpa(exp(-4), 10)命令可得的10位精确有效数字为2. 71828182&7. 389056099,54. 59815003,照可知，用级数法计算的有效数字分别为8, 4, 2, 0位。由此可以看岀，这个程序虽然原理上正确，但不好用。对不同的x,精度差别很大。英 他存在的问题有：这个程序不能用于X的元素群运算：当X为负数时，它成为交错级数，收敛很慢；此程序要做2次乘法，n很大时，乘法次数太多，计算速度很低；对不同的

7、x,要取不同的n 才能达到精度要求，因此n不应由用户输入，应该由软件按精度要求来选。正对上面的四个问题，可以采用下而四种方法改进：(1) 允许数组输入，改进输岀显示x二input (x二);n=input C n=,); y=ones(size (x);%输入原始数据，初始化 yfor i=l:ny=y+x. 4/prod(l:i); %循环相加sl=sprintf (* %13. Of, i); s2=sprintf %15. 8f, y); %将结果变为字符串disp(sl, s2)%显示end,执行此程序，输入x二12 4-4, n=10,结果为-3.0000000012.0000000

8、03.000000005.000000005.0000000022.500000005. 0000000013.00000000-5.6666666732.666666676.3333333323.666666675.0000000042.708333337.0000000034.33333333-3.5333333352.716666677.2666666742.866666672.1555555662.718055567.3555555648.55555556-1.0952381072.718253977.3809523851.806349210.5301587382.718278777.

9、3873015953.4317460392.718281537.3887125254.15414462-0.19223986102.718281807.3889947154.443104060.09671958(2) 可以利用exp(-x)=l/exp(x)来避免交错级数的计算：(3) 为了减少乘法次数，设一个中间变量z,它的初始值为z二ones (size (),把循环 体中的计算与句改为y=y+z; z=x*z/i;这样，求得的Z就是z二于是每个循环只需做一次乘法，计算整个级数只需n次乘 法。按这种计算，y的初始值改为y二zeros (size (x)(4) 为了按精度选择循环次数，不该使

10、用for循环，而用wh订e语句，它可以设置循环的 条件语句，通常可用y+z-ytol, tol是规左的允许误差只要相邻的两次y值之差大于tol,循环就继续进行，直到小于tol为I匕exp(x) = (exp(-)A当x较大时,exp(x)仍能很快收敛，还可以利用关系式k ,令xi=x/k. k通常取大于x而最接近x的2的幕，例如x二100,就取k二12&可以保证xl的绝对值小于1,这 时级数收敛得很快从练习中可以看岀,n取10时(即级数取10项)就能保证7位有效数，而 eXp(xl)，28可以化成* =(CXp( X1)2)2 )2,即exp(xl)的7次自乘，总共用17次乘法就 可完成两(1

11、00) = ( ( 100/128)2)2尸的计算,这既保证了精度,又提高了速度.例3编写任意函数展开为各阶泰勒级数的程序，并显示苴误差曲线.对于任意函数y=f(x),其泰勒展开式为/(x) = /() +- “) +(x “)2 + + 11 (x-” + Rn .2!n其中R3 为余项，也就是泰勒展开式的误差.MATLAB语句为fxs=input (输入y=f (x)的衣达式， s);%输入原始条件,fxs是字符串K=input (输入泰勒级数展开式的阶K);a=input (展开的位置a=);b=input (展开的区间半宽度b=，);x=linspace (a-b, a+b);%构成自

12、变虽数组，确定其长度和步长lx=length(x); dx=2*b/(lx-l);y=eval (fxs); %求岀y的准确值subplot (1, 2,1), plot (x, y,)，hold on %y 的准确值用点线绘出%求出a点的-阶导数，注意求导后数组长度减少1Dy=diff (y)/dx; Dya (1) =Dy (round (lx-1)/2);yt (1, :) =y (round (lx/2) +Dya(l)*(x-a); %求 y 的阶泰勒展开，绘图plot (x, yt (1,:)for k=2:KDy=diff(y, k)/(dxk); Dya(k) =Dy(roun

13、d(lx-k)/2); 滋求 3 点 k 阶导数yt(k, :)=yt(k-l, :) +Dya(k)/prod(1:k)*(x-a). k;弔求 y 的 k 阶导数 plot(x,yt(k, :);%绘图e(k, :)=y-yt(k, :);%求出 yt 的误差endtitle(fxs,的各阶泰初级数曲线),%注意如何组成标注的字符串grid, hold off, subplot (1, 2, 2)for k=l:K plot(x, e(k, :), hold on, end %绘制误差曲线title(fxs,的各阶泰勒级数误差曲线),grid, hold off执行此程序，输入fxs二co

14、s(x), K二5, a二0. 5, b二2,所得曲线见图3. 2(又变为误差曲线).读 者可以改变其坐标系范围以仔细观测最关心的部分，也可输入英他函数做验算，注意输入函 数应符合元素群运算规则.cos(x)的自阶事勒级数曲线cos(x)的各阶奉勒级数俣差曲线图31 y = COSX的泰勒级数及误差曲线严 1的值，可用symsum命令,相应的MATLAB代码为:工丄 计算级数心川2clear; syms k;simple(symsum(l/k2, 1, Inf) %simple 求解最简形式，Inf 为无穷人 结果为：ans 二l/6*pi2类似地可验证00 注:可用公式go川代码为：digi

15、ts(20);a=1.0; kk=l. 0;新殳置今后数值计算以20位相对精度进行%赋初值x i2A:为亠=J, = 12， 可以猜想有g汽5其中叫是正整数，请验证.来计算的近似值。如果要稱确到小数点后15位,相应的MATLABfor n=l:20, kk=kk/n; t a=a-kk;, endvpa(a, 17)%以17位相对精度给出a的值 结果为 w a 2.718281828459045.例5 (调和级数)1丄丄，丄自然数的倒数组成的数列2 3h 称为调和数列，由调JL 1JL 1和数列构成的级数心,2称为调和级数，我们把它的前项部分和k记为H（“）计算12JOO时H(n)和lnn,l

16、n(H+!)的值，并计算它们的差C(n) = H(n)-nnt c(n) = H(n)-ln(n +1),相应的曲佃代码为：H(1)=1; C(l)=l;for n=2:100, H(n) =H (n-l)+l/n;, %for 为循环语句c (n) =H (n) -log (n+1);, C (n) =H (n) -log (n);, end注意观测C(“)单调递减、c(“)单调递增，二者相互接近的现象。计算n = 10m(rn = 1,2,6)时C(n)和c(m)的值,注意观测C5)单调递减、旳)单 调递增，二者趋于同一极限的现象。并求出这个常数。极限C = lim(l + + - + + 丄 一 In