[time_fun.m]
function y=time_fun(t,T) % %
y=zeros(size(t));ii=find(t>=0.5 & t<=1.5);
y(ii)=ones(size(ii)).*(t(ii)-0.5);y(y==1.0)=0.5;
（2）
[t,y,S,aquad,bquad]=fzzyquad(0,2,15); % <1> AA=abs(A_sym);AA(AA<1e-10)=NaN; BB=abs(B_sym);BB(BB<1e-10)=NaN; A_quad=aquad(1:7) B_quad=bquad(1:7)
asq=abs(A_quad-A_sym)./AA bsq=abs(B_quad-B_sym)./BB A_quad =
0.2500 -0.3183 -0.0000 0.1061 0.0000 -0.0637 -0.0000 B_quad =
0 -0.2026 0.1592 0.0225 -0.0796 -0.0081 0.0531 asq =
1.0e-005 *
0 0.0105 NaN 0.0944 NaN 0.5316 NaN bsq =
1.0e-005 *
NaN 0.4866 0.1281 0.2949 0.3402 0.8191 0.7655
（3）图形显示截断余项后三角展开近似波形
SS=[S;y]';ribbon(t,SS);
xlabel('谐波次数\\rightarrow'),ylabel('t\\rightarrow')
view([46,38]),colormap(jet),shading flat,light,lighting gouraud
图 4.13-1
【例4.13.2.3-3】运用FFT，按式(4.13.2.2-2)和(4.13.2.2-3)，求上例时间函数的Fourier级数展开系数。 （1） [fzzyfft.m]
function [A,B,C,fn,t,w]=fzzyfft(T,M,Nf) %
31
if (nargin<2 | isempty(M));M=8;end if nargin<3;Nf=6;end N=2^M; % f=1/T; % w0=2*pi*f;
dt=T/N; % n=0:1:(N-1); % t=n*dt; %
w=time_fun(t,T); %
W=fft(w); % cn=W/N; % %
z_cn=find(abs(cn)<1.0e-10);
cn(z_cn)=zeros(length(z_cn),1); cn_SH=fftshift(cn); C=[cn_SH cn_SH(1)]; % A(1)=C(N/2+1);
A(2:N/2+1)=2*real(C((N/2+2):end)); B(2:N/2+1)=-2*imag(C((N/2+2):end));
if Nf>N/2;error(['第三输入宗量 Nf 应小于 ' int2str(N/2-1)]);end A(Nf+2:end)=[]; B(Nf+2:end)=[]; n1=-N/2:1:N/2; fn=n1*f;
（2）运行以下指令
% % <9> <10>
% <17>
% %
<23>
[A_fft,B_fft]=fzzyfft(2);
AA=abs(A_sym);AA(AA<1e-10)=NaN;BB=abs(B_sym);BB(BB<1e-10)=NaN; ast=abs(A_sym-A_fft)./AA bst=abs(B_sym-B_fft)./BB ast =
0 0.0001 NaN 0.0005 NaN 0.0013 NaN bst =
NaN 0.0001 0.0002 0.0005 0.0008 0.0013 0.0018
4.13.3 利用DFT计算一般连续函数的Fourier变换CFT 4.13.3.1 CFT与DFT之间的数学联系 4.13.3.2 MATLAB算法实现
【例4.13.3.2-1】运用FFT求取矩形脉冲w(t)??引起的混迭现象。 （1）
[cftbyfft.m]
function [AW,f]=cftbyfft(wt,t,flag) ?tbyfft.m
if nargin==2;flag=1;end
?1?00?t?1的谱，说明采样频率低else 32
N=length(t); % T=t(length(t))-t(1); % dt=T/N; %
W0=fft(wt); % <16> W=dt*W0; % df=1/T; % n=0:1:(N-1); %
if flag==0
n=-N/2:(N/2-1); W=fftshift(W); % end
f=n*df; % AW=abs(W); % if nargout==0
plot(f,AW);grid,xlabel('频率f');ylabel('|w(f)|') end
（2）
M=5; % <1> tend=1; %
T=10; % <3> N=2^M; % dt=T/N; % n=0:N-1; % t=n*dt; %
w=zeros(size(t,2),1);
Tow=find((tend-t)>0); % w(Tow,1)=ones(length(Tow),1); %
plot(t,w,'b','LineWidth',2.5),title('Time Waveform');xlabel('t --- >')
Time Waveform10.90.80.70.60.50.40.30.20.10024t --- >6810图4.13-2
[AW,f]=cftbyfft(w,t,0); ff=f+eps;
AWW=abs(sin(pi*ff)./(pi*ff)); plot(f,AW,'b-',ff,AWW,'r:')
title('Aliasing caused by undersampling')
xlabel('f --- >');ylabel('|W(f)|'),legend('by FFT','Theoretical')
33
Aliasing caused by undersampling1.41.210.80.60.40.20-2by FFT Theoretical|W(f)|-1.5-1-0.50f --- >0.511.52图4.13-3
4.14 常微分方程
4.14.1 初值常微分方程的解算指令 4.14.1.1 解ODE的基本机理 4.14.1.2 solver解算指令的使用说明 4.14.2 ODE解算指令的使用演示 4.14.2.1 解算指令简洁格式使用示例
【例4.14.2.1-1】采用最简洁格式的ODE文件和解算指令，研究围绕地球旋转的卫星轨道。 （1）
（2）
（3） [dYdt.m]
function Yd=DYdt(t,Y) % t % Y
global G ME % xy=Y(1:2);Vxy=Y(3:4); % r=sqrt(sum(xy.^2));
Yd=[Vxy;-G*ME*xy/r^3]; %
（4）
global G ME % <1>
G=6.672e-11;ME=5.97e24;vy0=4000;x0=-4.2e7;t0=0;tf=60*60*24*9; tspan=[t0,tf]; % Y0=[x0;0;0;vy0]; %
[t,YY]=ode45('DYDt',tspan,Y0);% <8> X=YY(:,1); % Y=YY(:,2); %
plot(X,Y,'b','Linewidth',2); hold on
34
axis('image') [XE,YE,ZE] = sphere(10); RE=0.64e7; XE=RE*XE;YE=RE*YE;ZE=0*ZE; mesh(XE,YE,ZE),hold off x 1086420-2-4-6-807%
% % % % 5101520x 107图 4.14-2
【例4.14.2.1-2】上例中，程序间的参数（如G和ME）传送，是依靠全局变量形式实现的。一般说来，编写程序时，应尽量少用全局变量，以免引起混乱。本例演示参数如何在指令间直接传送。 （1）
[DYDt2.m]
function Yd=DYDt2(t,Y,flag,G,ME) % flag %
switch flag
case '' %
X=Y(1:2);V=Y(3:4);r=sqrt(sum(X.^2));Yd=[V;-G*ME*X/r^3]; otherwise
error(['Unknown flag ''' flag '''.']); end
（2）
[t,YY]=ode45('DYDt2',tspan,Y0,[],G,ME);
4.14.2.2 解算指令较复杂格式的使用示例
【例4.14.2.2-1】带事件设置的ODE文件及主程序编写演示。本例将以较高精度计算卫星经过近地点和远地点的时间，并在图上标志。 （1）
[DYDt3.m]
function varargout=DYDt3(t,Y,flag,G,ME,tspan,Y0) % % % %
35