数值分析第三章上机实习题

第一题

没什么好说的.用polyfit拟合就可以了.

clear;f=@(x) 1./(1+25*x.*x);x = -1 + 0.2*(0:10);y = f(x);P = polyfit(x,y,3);u = -1:0.01:1;v_1 = f(u);v_2 = polyval(P,u);plot(x,y,'ro',u,v_1,'k-',u,v_2,'b--');

第二题

先做3次4次的拟合.

clear;x = [0,.1,.2,.3,.5,.8,1];y = [1,.41,.5,.61,.91,2.02,2.46];s_3 = polyfit(x,y,3);s_4 = polyfit(x,y,4);y_3 = polyval(s_3,x);y_4 = polyval(s_4,x);u = 0:0.01:1;v_3 = polyval(s_3,u);v_4 = polyval(s_4,u);plot(u,v_3,'r--',u,v_4,'b--',x,y,'ko',x,y_3,'ro',x,y_4,'bo');

其实拟合的已经很不错了= =,但是题目要求再用一种曲线拟合一次,那就先对数据求对数,发现大致在一条直线左右,然后进行一次拟合.

lny = log(y);s_0 = polyfit(x,lny,1);vv = polyval(s_0,u);plot(x,lny,'ko',u,vv,'b-');

最后再把数据转换回来.

y_0 = exp(polyval(s_0,x));v_0 = exp(polyval(s_0,u));plot(x,y,'ko',x,y_0,'bo',u,v_0,'b--');

第三题

使用快速傅里叶变换确定函数f(x)=x^2*cosx在[-π,π]上的16次三角插值多项式

matlab中的fft和ifft和书中的定义是有差别的,可以先输入help fft看一下公式

先拿书上91页的例13讲一下

如果是按照书上定义的ck,那么应该用ifft函数

clear;syms x;fx = x^4 - 3*x^3 + 2*x^2 - tan(x*(x-2));u = linspace(0,2,200);v = subs(fx,u);j = 0:7;fi = subs(fx,j/4);c = 8*ifft(fi);an = real(c.*exp(-1i*pi*j))/4;bn = imag(c.*exp(-1i*pi*j))/4;s = an(1)/2 + an(5)*cos(4*x);for ix = 2:4    s = s + an(ix)*cos((ix-1)*x) + bn(ix)*sin((ix-1)*x);endv0 = subs(s,pi*(u-1));plot(u,v,'k-',u,v0,'b-');

如果要用fft函数,那需要小小的改一下c,an,bn,如下:

c = fft(fi);an = real(c.*exp(-1i*pi*j))/4;bn = -imag(c.*exp(-1i*pi*j))/4;

回到这道题,如果用ifft

clear;syms x;fx = x*x*cos(x);j = 0:31;u = linspace(-pi,pi,200);v = subs(fx,u);xx = -pi + pi*j/16;yy = subs(fx,xx);c = 32*ifft(yy);an = real(c.*exp(-1i*pi*j))/16;bn = imag(c.*exp(-1i*pi*j))/16;s = an(1)/2 + an(17)*cos(16*x);for ix = 2:16    s = s + an(ix)*cos((ix-1)*x) + bn(ix)*sin((ix-1)*x);endv0 = subs(s,u);plot(u,v0,'r-',u,v,'k-');

fft,改c,an,bn

c = fft(yy);an = real(c.*exp(-1i*pi*j))/16;bn = -1*imag(c.*exp(-1i*pi*j))/16;

在网上还搜到了这样的一个函数

function [an,bn,f]=fseries(fx,x,n,a,b)if nargin==3    a=-pi;    b=pi;endl=(b-a)/2;if a+b    fx=subs(fx,x,x+l+a);endan=int(fx,x,-l,l)/l;bn=[];f=an/2;for ii=1:n    ann=int(fx*cos(ii*pi*x/l),x,-l,l)/l;    bnn=int(fx*sin(ii*pi*x/l),x,-l,l)/l;    an=[an,ann];    bn=[bn,bnn];    f=f+ann*cos(ii*pi*x/l)+bnn*sin(ii*pi*x/l);endif a+b    f=subs(f,x,x-l-a);end

这个根本就不是离散傅里叶变换,而是傅里叶变换,其中的int为求积分,所以这个函数运行速度很慢!