基于MATLAB的幂级数求和与展开(Taylor和Fourier算法) | 您所在的位置:网站首页 › matlab设计一个函数求和 › 基于MATLAB的幂级数求和与展开(Taylor和Fourier算法) |
目录 一. Taylor幂级数展开 1.1 单变量 例题1 例题2 1.2 多变量 例题3 二. Fourier级数展开 例题4 例题5 三. 级数求和 例题6 例题7 例题8 例题9 一. Taylor幂级数展开 1.1 单变量在x=0点附近的Taylor幂级数如下: 在x=a点展开: MATLAB格式: %x=0进行Taylor幂级数展开 taylor(f,x,k) %x=a进行Taylor幂级数展开 taylor(f,x,k,a) 例题1求函数Taylor幂级数展开的前9项;并求关于x=2和x=a的Taylor幂级数展开。 解: MATLAB代码: clc;clear; syms x a; f=sin(x)/(x^2+4*x+3); y1=taylor(f,x,'order',9); pretty(y1) taylor(f,x,2,'order',9) taylor(f,x,a,'order',5)运行结果: ans =sin(2)/15 - ((131623*cos(2))/35880468750 + (875225059*sin(2))/34445250000000)*(x - 2)^8 + (x - 2)*(cos(2)/15 - (8*sin(2))/225) - ((623*cos(2))/11390625 + (585671*sin(2))/2733750000)*(x - 2)^6 + ((262453*cos(2))/19136250000 + (397361*sin(2))/5125781250)*(x - 2)^7 - (x - 2)^2*((8*cos(2))/225 + (127*sin(2))/6750) + (x - 2)^3*((23*cos(2))/6750 + (628*sin(2))/50625) + (x - 2)^4*((28*cos(2))/50625 - (15697*sin(2))/6075000) + (x - 2)^5*((203*cos(2))/6075000 + (6277*sin(2))/11390625) ans =sin(a)/(a^2 + 4*a + 3) - (cos(a)/(a^2 + 4*a + 3) - (sin(a)*(2*a + 4))/(a^2 + 4*a + 3)^2)*(a - x) - (a - x)^2*(sin(a)/(2*(a^2 + 4*a + 3)) + (cos(a)*(2*a + 4))/(a^2 + 4*a + 3)^2 - (sin(a)*((2*a + 4)^2/(a^2 + 4*a + 3)^2 - 1/(a^2 + 4*a + 3)))/(a^2 + 4*a + 3)) - (a - x)^3*((sin(a)*((2*a + 4)/(a^2 + 4*a + 3)^2 - (((2*a + 4)^2/(a^2 + 4*a + 3)^2 - 1/(a^2 + 4*a + 3))*(2*a + 4))/(a^2 + 4*a + 3)))/(a^2 + 4*a + 3) - cos(a)/(6*(a^2 + 4*a + 3)) + (cos(a)*((2*a + 4)^2/(a^2 + 4*a + 3)^2 - 1/(a^2 + 4*a + 3)))/(a^2 + 4*a + 3) + (sin(a)*(2*a + 4))/(2*(a^2 + 4*a + 3)^2)) + (a - x)^4*(sin(a)/(24*(a^2 + 4*a + 3)) - (sin(a)*(((2*a + 4)^2/(a^2 + 4*a + 3)^2 - 1/(a^2 + 4*a + 3))/(a^2 + 4*a + 3) + ((2*a + 4)*((2*a + 4)/(a^2 + 4*a + 3)^2 - (((2*a + 4)^2/(a^2 + 4*a + 3)^2 - 1/(a^2 + 4*a + 3))*(2*a + 4))/(a^2 + 4*a + 3)))/(a^2 + 4*a + 3)))/(a^2 + 4*a + 3) + (cos(a)*((2*a + 4)/(a^2 + 4*a + 3)^2 - (((2*a + 4)^2/(a^2 + 4*a + 3)^2 - 1/(a^2 + 4*a + 3))*(2*a + 4))/(a^2 + 4*a + 3)))/(a^2 + 4*a + 3) + (cos(a)*(2*a + 4))/(6*(a^2 + 4*a + 3)^2) - (sin(a)*((2*a + 4)^2/(a^2 + 4*a + 3)^2 - 1/(a^2 + 4*a + 3)))/(2*(a^2 + 4*a + 3))) 例题2对sinx进行Taylor幂级数展开,并观察不同阶次的近似效果。 解: MATLAB代码: clc;clear; x0=-2*pi:0.01:2*pi; y0=sin(x0); syms x; y=sin(x); plot(x0,y0,'r-'),axis([-2*pi,2*pi,-1.5,1.5]); hold on for n=[8:2:16] p=taylor(y,x,'order',n), y1=subs(p,x,x0); line(x0,y1); end运行结果: p =- x^7/5040 + x^5/120 - x^3/6 + x p =x^9/362880 - x^7/5040 + x^5/120 - x^3/6 + x p =- x^11/39916800 + x^9/362880 - x^7/5040 + x^5/120 - x^3/6 + x p =x^13/6227020800 - x^11/39916800 + x^9/362880 - x^7/5040 + x^5/120 - x^3/6 + x p =- x^15/1307674368000 + x^13/6227020800 - x^11/39916800 + x^9/362880 - x^7/5040 + x^5/120 - x^3/6 + x 1.2 多变量多变量函数在的Taylor幂级数展开如下: 例题3求z=f(x,y)的各种Taylor幂级数展开。 解: MATLAB代码: clc;clear; syms x y a; f=(x^2-2*x)*exp(-x^2-y^2-x*y); F1=taylor(f,[x,y],[0,0],'order',8) F2=taylor(f,[x,y],[1,a],'order',3)运行结果: F1 = x^7/3 + x^6*y + x^6/2 + 2*x^5*y^2 + x^5*y - x^5 + (7*x^4*y^3)/3 + (3*x^4*y^2)/2 - 2*x^4*y - x^4 + 2*x^3*y^4 + x^3*y^3 - 3*x^3*y^2 - x^3*y + 2*x^3 + x^2*y^5 + (x^2*y^4)/2 - 2*x^2*y^3 - x^2*y^2 + 2*x^2*y + x^2 + (x*y^6)/3 - x*y^4 + 2*x*y^2 - 2*x F2 = (exp(- a^2 - a - 1) - exp(- a^2 - a - 1)*((a/2 + 1)*(a + 2) - 1))*(x - 1)^2 - exp(- a^2 - a - 1) - exp(- a^2 - a - 1)*(a - y)^2*((2*a + 1)*(a + 1/2) - 1) + exp(- a^2 - a - 1)*(a + 2)*(x - 1) - exp(- a^2 - a - 1)*(2*a + 1)*(a - y) + exp(- a^2 - a - 1)*(a - y)*(x - 1)*((2*a + 1)*(a/2 + 1) + (a + 2)*(a + 1/2) - 1) 二. Fourier级数展开给定f(x)以及如下条件: Fourier级数,如下: 编写Fourier级数的解析函数: function [A,B,F]=fseries(f,x,n,a,b) if nargin==3,a=-pi; b=pi; end L=(b-a)/2; if a+b, f=subs(f,x,x-L-a); end %变量区域互换 A=int(f,x,-L,L)/L; B=[]; F=A/2; %计算a0 for i=1:n an=int(f*cos(i*pi*x/L),x,-L,L)/L; bn=int(f*cos(i*pi*x/L),x,-L,L)/L; A=[A,an]; B=[B,bn]; F=F+an*cos(i*pi*x/L)+bn*sin(i*pi*x/L); end if a+b, F=subs(F,x,x+L+a); end %换回变量区域 end 例题4求y的Fourier级数展开。 解: MATLAB代码: clc;clear; syms x; f=x*(x-pi)*(x-2*pi); [A,B,F]=fseries(f,x,6,0,2*pi) function [A,B,F]=fseries(f,x,n,a,b) if nargin==3,a=-pi; b=pi; end L=(b-a)/2; if a+b, f=subs(f,x,x-L-a); end %变量区域互换 A=int(f,x,-L,L)/L; B=[]; F=A/2; %计算a0 for i=1:n an=int(f*cos(i*pi*x/L),x,-L,L)/L; bn=int(f*sin(i*pi*x/L),x,-L,L)/L; A=[A,an]; B=[B,bn]; F=F+an*cos(i*pi*x/L)+bn*sin(i*pi*x/L); end if a+b, F=subs(F,x,x+L+a); end %换回变量区域 end运行结果: A =[ -16*pi^3, 24*pi, -6*pi, (8*pi)/3, -(3*pi)/2, (24*pi)/25, -(2*pi)/3] B = [ -(12*pi - 24*pi^3)/pi, ((3*pi)/2 - 12*pi^3)/pi, -((4*pi)/9 - 8*pi^3)/pi, ((3*pi)/16 - 6*pi^3)/pi, (24*pi^2)/5 - 12/125, (pi/18 - 4*pi^3)/pi] F = (sin(x)*(12*pi - 24*pi^3))/pi - 24*pi*cos(x) - 8*pi^3 - 6*pi*cos(2*x) - (8*pi*cos(3*x))/3 - (3*pi*cos(4*x))/2 - (24*pi*cos(5*x))/25 - (2*pi*cos(6*x))/3 - sin(5*x)*((24*pi^2)/5 - 12/125) + (sin(2*x)*((3*pi)/2 - 12*pi^3))/pi + (sin(3*x)*((4*pi)/9 - 8*pi^3))/pi + (sin(4*x)*((3*pi)/16 - 6*pi^3))/pi + (sin(6*x)*(pi/18 - 4*pi^3))/pi 例题5对区间内的方波信号,当时,y=1,否则y=-1。对它进行Fourier级数拟合,并观察用多少项能有较好的拟合效果。 解: MATLAB代码: clc;clear; syms x; f=abs(x)/x; %定义方波信号 xx=[-pi:pi/200:pi]; xx=xx(xx~=0); xx=sort([xx,-eps,eps]); %剔除零点 yy=subs(f,x,xx); plot(xx,yy,'r-.'), %绘制出理论值 hold on %保持坐标系 for n=2:20 [a,b,f1]=fseries(f,x,n); y1=subs(f1,x,xx); plot(xx,y1) end function [A,B,F]=fseries(f,x,n,a,b) if nargin==3,a=-pi; b=pi; end L=(b-a)/2; if a+b, f=subs(f,x,x-L-a); end %变量区域互换 A=int(f,x,-L,L)/L; B=[]; F=A/2; %计算a0 for i=1:n an=int(f*cos(i*pi*x/L),x,-L,L)/L; bn=int(f*sin(i*pi*x/L),x,-L,L)/L; A=[A,an]; B=[B,bn]; F=F+an*cos(i*pi*x/L)+bn*sin(i*pi*x/L); end if a+b, F=subs(F,x,x+L+a); end %换回变量区域 end运行结果: 三. 级数求和给定: 符号工具箱中有提供: S=symsum(fk,k,k0,kn) 例题6计算求和: 解: MATLAB代码: clc;clear; format long; %数值计算 sum(2.^[0:63]) %把2定义为符号量,计算会更准确 sum(sym(2).^[0:63]) %工具箱 syms k; symsum(2^k,0,63)运行结果: ans =1.844674407370955e+19 ans =18446744073709551615 ans =18446744073709551615 例题7试求无穷级数的和 解: MATLAB代码: clc;clear; %符号工具箱 syms n; s=symsum(1/((3*n-2)*(3*n+1)),n,1,inf) %数值计算 format long %以长型方式显示得出的结果 m=1:1000000; s1=sum(1./((3*m-2).*(3*m+1))) %双精度有效位16运行结果: s =1/3 s1 =0.333333222222265 例题8求解: 解: MATLAB代码: clc;clear; syms n x; s1=symsum(2/((2*n+1)*(2*x+1)^(2*n+1)),n,0,inf)运行结果: s1 =piecewise(1 < abs(2*x + 1), 2*atanh(1/(2*x + 1))) 例题9求: 解: MATLAB代码: clc;clear; syms m n; limit(symsum(1/m,m,1,n)-log(n),n,inf); vpa(ans,70) %显示70位有效数字运行结果: ans = 0.5772156649015328606065120900824024310421593359399235988057672348848677 |
CopyRight 2018-2019 实验室设备网 版权所有 |