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

(1 x) 1 x ( 1) x 2 , x (1,1). 2!
作图观测其展开式的前几项多项式逼近原函数的情况。
例 2 利用幂级数计算指数函数。指数函数可展开为幂级数
e x 1 x x 2 x3 x n , x (,),
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
数
3!
3! 5! 的 图 形 ， 观 测 这 些 多 项 式 函 数 的 图 形 向
y sin x 的图形的逼近的情况。例如，在区间 [0, ] 上作函数 y sin x 与多项式函数
y x, y x x3 , y x x3 x5
3!
3! 5! 图形的 MATLAB 代码为：
(4) 为了按精度选择循环次数,不该使用 for 循环,而用 while 语句,它可以设置循环的条件
3
语句,通常可用 y+z-y>tol,tol 是规定的允许误差.只要相邻的两次 y 值之差大于 tol,循环就继续 进行,直到小于 tol 为止.
exp( x) (exp( x ))k
当 x 较大时,exp(x)仍能很快收敛,还可以利用关系式
的值,可用
symsum
命令,相应的
MATLAB
代码为：
>>clear; syms k;
>>simple(symsum(1/k^2,1,Inf)) 结果为:
ans =1/6*pi^2 类似地可验证
%simple 求解最简形式,Inf 为无穷大
1
n4
n1
4 , 90
1 n6
n1

6 .
2! 3!
n!
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
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
>>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(又变为误差曲线).读者可 以改变其坐标系范围以仔细观测最关心的部分,也可输入其他函数做验算,注意输入函数应符 合元素群运算规则.
执行此程序，分别带入 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, e2 , e4 , e4 的 10
>>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 sin x 的图形。
1.5
1
0.5
0
0
1

f (n) (a) (x a)n n!
Rn (x).
其中
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
完成 exp(100) (((exp(100 /128))2 )2 )2 的计算,这既保证了精度,又提高了速度.
例 3 编写任意函数展开为各阶泰勒级数的程序,并显示其误差曲线.对于任意函数 y=f(x),其泰勒展开式为
f (x)
f (a)
f '(a)(x a)
f
(2) (a) (x a)2 2!
（3）为了减少乘法次数，设一个中间变量 z，它的初始值为 z=ones(size(x)),把循环体中
的计算与句改为
y=y+z; z=x.*z/i;
这样，求得的 z 就是 z=x.^i/i!,于是每个循环只需做一次乘法，计算整个级数只需 n 次乘法。
按这种计算，y 的初始值改为 y=zeros(size(x))
945
1
可以猜想有 n1 n 2k
2k , k 1,2,,
mk
其中 m k 是正整数,请验证.


1

e
注：可用公式 n0 n! 来计算 e 的近似值。如果要精确到小数点后 15 位，相应的
k ,令 x1=x/k.k 通
常取大于 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 次乘法就可
位精确有效数字为 2.718281828, 7.389056099, 54.59815003, .1831563889e-1
对照可知，用级数法计算的有效数字分别为 8，4，2，0 位。 由此可以看出，这个程序虽然原理上正确，但不好用。对不同的 x,精度差别很大。其
他存在的问题有： 这个程序不能用于 x 的元素群运算；当 x 为负数时，它成为交错级数，收敛很慢；此
4
cos(x)的 各 阶 泰 勒 级 数 曲 线 2
1.5
1
0.5
0
-0.5
-1
-1.5
-2
-2
0
2
4
cos(x)的 各 阶 泰 勒 级 数 误 差 曲 线 1.2
1
0.8
0.6
0.4
0.2
0
-0.2
-2
0
2
4
图 3.1 y cos x 的泰勒级数及误差曲线
例4
计算级数
n1
1 n2
.
可以用 help symsum, help taylor 查阅有关这些命令的详细信息
例 1 先用 taylor 命令观测函数 y sin x 的 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)