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

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

(7页)

本资源提供全文预览，点击全文预览即可全文预览,如果喜欢文档就下载吧，查找使用更方便哦！

14.90 积分

实验七　matlab求解级数有关计算    1.级数的基本概念常数项级数：称用加号将数列的项连成的式子为(常数项)无穷级数，简记为。称级数前项构成的和为级数的部分和。若，则称级数收敛，其和为。Taylor级数：设函数在包含的区域内具有各阶导数，则称幂级数为函数在的Taylor级数，当时称为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命令观测函数的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 =xans =xans =x-1/6*x^3ans =x-1/6*x^3ans =x-1/6*x^3+1/120*x^5然后在同一坐标系里作出函数和它的Taylor展开式的前几项构成的多项式函数的图形，观测这些多项式函数的图形向的图形的逼近的情况。例如，在区间上作函数与多项式函数图形的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，其中实线表示函数的图形。图3.1 的泰勒级数类似地，根据函数的Taylor级数作图观测其展开式的前几项多项式逼近原函数的情况。例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)命令可得的10位精确有效数字为2.718281828,  7.389056099,  54.59815003,  .1831563889e-1对照可知，用级数法计算的有效数字分别为8，4，2，0位。由此可以看出，这个程序虽然原理上正确，但不好用。对不同的x,精度差别很大。其他存在的问题有：这个程序不能用于x的元素群运算；当x为负数时，它成为交错级数，收敛很慢；此程序要做次乘法，n很大时，乘法次数太多，计算速度很低；对不同的x,要取不同的n才能达到精度要求，因此n不应由用户输入，应该由软件按精度要求来选。正对上面的四个问题，可以采用下面四种方法改进：(1)允许数组输入，改进输出显示x=input('x='); n=input('n='); y=ones(size(x));  %输入原始数据，初始化yfor 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语句,它可以设置循环的条件语句,通常可用y+z-y>tol,tol是规定的允许误差.只要相邻的两次y值之差大于tol,循环就继续进行,直到小于tol为止.当x较大时,exp(x)仍能很快收敛,还可以利用关系式,令x1=x/k.k通常取大于x而最接近x的2的幂,例如x=100,就取k=128,可以保证x1的绝对值小于1,这时级数收敛得很快..从练习中可以看出,n取10时(即级数取10项)就能保证7位有效数,而可以化成,即exp(x1)的7次自乘,总共用17次乘法就可完成的计算,这既保证了精度,又提高了速度.例3       编写任意函数展开为各阶泰勒级数的程序,并显示其误差曲线.对于任意函数y=f(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(又变为误差曲线).读者可以改变其坐标系范围以仔细观测最关心的部分,也可输入其他函数做验算,注意输入函数应符合元素群运算规则.图3.1 的泰勒级数及误差曲线例4       计算级数的值,可用symsum命令,相应的MATLAB代码为：>>clear; syms k;>>simple(symsum(1/k^2,1,Inf))  %simple求解最简形式,Inf为无穷大结果为:    ans =1/6*pi^2类似地可验证可以猜想有其中是正整数,请验证.注：可用公式来计算的近似值。如果要精确到小数点后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的值结果为例5 (调和级数 )  自然数的倒数组成的数列称为调和数列，由调和数列构成的级数称为调和级数，我们把它的前项部分和记为。计算时和的值，并计算它们的差，，相应的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注意观测单调递减、单调递增，二者相互接近的现象。  计算时和的值，注意观测单调递减、单调递增，二者趋于同一极限的现象。并求出这个常数。  极限称为欧拉(Euler)常数，可以证明它是一个无理数。显然，而由于当时趋于0，故趋于同一个极限。 习题16-71. 用taylor命令观测函数的Maclaurin展开式的前几项, 然后在同一坐标系里作出函数和它的Taylor展开式的前几项构成的多项式函数的图形，观测这些多项式函数的图形向的图形的逼近的情况(1)   (2)  (3) (4)    (5)     (6) 2. 求公式中的数的值.3. 利用公式来计算的近似值。精确到小数点后100位,这时应计算到这个无穷级数的前多少项?请说明你的理由.4. 用练习3中所用观测法判断下列级数的敛散性(1)   (2)   (3)   (4) (5)   (6)   (7)   (8) 关 键 词： 实验 matlab 求解 级数 有关 计算

 天天文库所有资源均是用户自行上传分享，仅供网友学习交流，未经上传用户书面授权，请勿作他用。