傅里叶逼近

《MATLAB程序设计实践》

科学计算的算法及其举例应用：

(1)、傅里叶逼近：

1.用法说明：对于连续周期函数，只要计算出其傅里叶展开级数即可，在Matlab中编程实现的连续函数的傅里叶逼近法函数为：FZZ。

功能：用傅里叶级数逼近已知的连续周期函数。

调用格式：[A0，A，B]=FZZ(func,T,n).

其中，func为已知函数；

T为已知函数的周期；

N为展开级数的项数；

A0为展开后的常数项；

A为展开后的余弦项系数；

B为展开后的正弦项系数。

2.源程序代码：

function [A0,A,B]=FZZ(func,T,n)

syms t;

func=subs(sym(func),sym('x'),sym('t'));

A0=int(sym(func),t,-T/2,T/2)/T;

for(k=1:n)

A(k)=int(func*cos(2*pi*k*t/T),t,-T/2,T/2)*2/T;

A(k)=vpa(A(k),4);

B(k)= int(func*sin(2*pi*k*t/T),t,-T/2,T/2)*2/T;

B(k)= vpa(B(k),4);

end

3.举例说明：傅里叶逼近应用实例。用傅里叶级数(取5项)逼近函数x，输出系数值。

>> [A0,A,B]=FZZ('x',2*pi,5)

A0 =

0

A =

[ 0., 0., 0., 0., 0.]

B =

[ 2., -1., .6667, -.5000, .4000]

3.流程图(1)

二、科学计算和工程实际问题：

(1)设单自由度阻尼系统的质量M=1kg,弹簧刚度系数K=100N/m，速度阻尼系数c=4N*s/m,求它在如下外力下的强迫振动，得出t≦1.2s的波形。

f=t/0.015(0≦t≦0.15),

f=10(0.15≦t≦1.2)

系统的受力平衡表达式：

f-c*v-K*x=M*a; (1)

其中x为位移，v为速度，a为加速度

用dx/dt代替v，得：

f-4*dx/dt-100*x=1*d2x/dt2; (2)

运用函数分段求解位移的表达式，然后用plot画图

1.源程序：

>> dsolve('D2x=t/0.015-4*Dx-100*x','x(0)=0','Dx(0)=0')

ans =

-23/900*exp(-2*t)*sin(4*6^(1/2)*t)*6^(1/2)+2/75*exp(-2*t)*cos(4*6^(1/2)*t)-2/75+2/3*t

>> t1=0:0.01:0.15;

>> x1=-23/900*exp(-2.*t1).*sin(4*6^(1/2).*t1).*6^(1/2)+2/75*exp(-2.*t1).*cos(4*6^(1/2).*t1)-2/75+2/3.*t1;

>> subs(x1,t1,0.15)

ans =

Columns 1 through 11

0 0.0000 0.0001 0.0003 0.0007 0.0013 0.0022 0.0035 0.0051 0.0071 0.0096

Columns 12 through 16

0.0125 0.0160 0.0199 0.0243 0.0292

>> subs(diff(x1),t1,0.15)

ans =

Columns 1 through 11

0.0000 0.0001 0.0002 0.0004 0.0006 0.0009 0.0013 0.0016 0.0020 0.0025 0.0029

Columns 12 through 15

0.0034 0.0039 0.0044 0.0049

>> dsolve('D2x=10-4*Dx-100*x','x(0.15)=0.0292','Dx(0.15)=0.0049')

ans =

-1/240000*exp(-2*t