切比雪夫级数

1、算法说明：当一个连续函数定义在区间-1,1上时，它可以展开成切比雪夫级数。即：其中为次切比雪夫多项式，具体表达可通过递推得出： ， 它们之间满足如下的正交关系： 在实际应用中，可根据所需的精度来截取有限的项数，切比雪夫级数中的系数由下式决定： 在MATLAB中编程实现的切比雪夫逼近法函数为：Chebyshev。功能：用切比雪夫多项式逼近已知函数。调用格式：其中，y为已知函数； k为逼近已知函数所需项数； f是求得的切比雪夫逼近多项式在x0处的逼近值。程序源代码（m文件）：function f = Chebyshev(y,k,x0)%用切比雪夫多项式逼近已知函数%已知函数：y%逼近已知函数所需项

2、数：k%逼近点的x坐标：x0%求得的切比雪夫逼近多项式或在x0处的逼近值：f syms t;T(1:k+1) =t;T(1) = sym(1);T(2) = t;c(1:k+1) = sym(0); c(1)=int(subs(y,findsym(sym(y),sym(t)*T(1)/sqrt(1-t2),t,-1,1)/pi;c(2)=2*int(subs(y,findsym(sym(y),sym(t)*T(2)/sqrt(1-t2),t,-1,1)/pi;f = c(1)+c(2)*t; for i=3:k+1 T(i) = 2*t*T(i-1)-T(i-2); c(i) = 2*int(

3、subs(y,findsym(sym(y),sym(t)*T(i)/sqrt(1-t2),t,-1,1)/2; f = f + c(i)*T(i); f = vpa(f,6); if(i=k+1) if(nargin = 3) f = subs(f,t,x0); else f = vpa(f,6); end endend应用实例：切比雪夫应用实例。用切比雪夫公式（取6项）逼近函数，并求当x=0.5时的函数值。解：利用程序求解方程，在MATLAB命令窗口中输入： Chebyshev(1/(2-x),6) %调用创建的函数euler，输出切比雪夫多项式的6个项再在MATLAB命令窗口中输入： Ch

4、ebyshev(1/(2-x),6,0.5) %调用创建的函数euler，输出当x=0.5时的函数值输出结果：由第一项和第二项多项式及其系数计算一级切比雪夫逼近值根据前两项多项式计算其系数c(1)，c(2)定义sym型变量t以及切比雪夫多项式矩阵T(n)，并规定前两项为1，t开始流程图： i从3到k+1 是由T(i-1)和T(i-2)计算T(i) 否由T(i)计算其系数c(i)结束输出结果，并设置其为六级精度求出逼近式在x0处的值输入三项i=k+1由T(i)和c(i)计算i级切比雪夫逼近值。 是 是 否二、拉压杆系的静不定问题。由n根杆（CB1，CB2.CBn）组成的桁架结构如图2-1所示，受

5、力P作用，各杆的横截面积分别为Ai，材料弹性模量为E，长度为Li，求各杆的轴力Ni以及节点C处的位移。 图2-11、假设：由P的作用，C点移动到C点；C到C的水平距离为dx，垂直距离为dy；各杆与水平面夹角分别为1，2in；每根杆长度的变化为dL1，dL2dLidLn；力P与水平面夹角为。2、思路：（1）如图2-2所示，以CBi为例说明各杆的几何关系，由于，0，因此dLi =BiC-BiCBiC-BiD=DC，过C点作DC的平行线CG，因，故，所以CF=dxcosi，FG=dysini,故可得几何方程。图2-2 （2）根据图2-3所示各轴力以及外力x,y方向合力为0，建立平衡方程、：图2-3

6、（3）由公式、共有n+2个公式，求解n个轴力，以及A点位移dx、dy，建立如下所示的线性方程组：.（4）建立Pcos,Psin，0，0，00的常数矩阵，以及如下所示的系数矩阵（5）再用求逆法求解此线性方程组，即用常数矩阵除以系数矩阵，得出结果。3、源程序：（文件名称为main）clear;clc;Ei=input(请输入各杆的刚度:(注意用括起来); %输入刚度矩阵EiLi=input(请输入各杆的长度:(注意用括起来) );%输入杆的长度矩阵LiAi=input(请输入各杆的横截面积:(注意用括起来) );%输入杆的横截面积矩阵Aiai=input(请输入各杆与水平面的夹角:(注意用括起来)

7、 );%输入杆与水平面的夹角矩阵aiP=input(请输入外力P: );%输入外加力Pa=input(请输入P与水平面的夹角: )%输入外加力P与x的夹角n1=length(Ei);n2=length(Li);n3=length(Ai);if(n1=n2|n2=n3|n1=n3) disp(输入数据错误)else n=n1; end%判断数据大小是否一致Ki=Li./(Ei.*Ai);C=zeros(n+2,1);C(1,1)=P*cos(a);C(2,1)=P*sin(a);C(3:n+2,1)=zeros(n,1);%建立方程组等号右边常数的矩阵D=zeros(n+2,n+2);D(1,:

8、)=cos(ai),0,0;D(2,:)=sin(ai),0,0;for(i=1:n) D(i+2,i)=Ki(i);endD(3:n+2,n+1)=(-cos(ai);D(3:n+2,n+2)=(-sin(ai);%建立方程组系数矩阵x=DC;x=x;%求解该线性方程组，得出C点位移以及每根杆的轴力disp(节点在x、y方向上的位移分别：)x(n+1:n+2)disp(各杆的轴力分别为：)x(1:n)%输出结果4、流程图：开始输入所需数据：各杆的刚度Ei、横截面积Ai、长度Li，以及夹角i，外加力P以及其与x方向夹角输入Ei，Ai，Li，i维度是否一致是否建立系数矩阵D，以及常数矩阵C输出：

9、输入数据错误用求逆法解线性微分方程输出C点位移，以及各杆轴力结束5、程序举例应用：设三根杆组成的支架如图2-4所示，挂一重物P=3000N。设L=3m，各杆的横截面积分别为：A1=15010-6m2，A2=20010-6m2，A3=30010-6m2，材料的弹性模量均为E=200109N/m2，求各杆所受力的大小以及C点位移图2-5解：运行main文件（即前文中写的源程序）；输入题中所给数据； 200e9,200e9,200e9; 3/sin(pi/3),3/sin(pi/2),3/sin(pi/4); 150e-6,200e-6,300e-6; pi/3,pi/2,3*pi/4;3000;0

10、;得出结果；（1）；一、流程图：建立函数文件输入MATLAB指令 若不成立 修改结果二、源程序代码：建立被积函数文件funofuno.mfunction f=funo(x)f=(1./(2*pi).0.5).*exp(-x.2./2);计算积分，在MATLAB命令窗中输入： quad(funo,0,1)三、结果：ans =0.3413（2）；一、流程图：建立函数文件输入MATLAB指令 若不成立 修改结果二、源程序代码：计算积分，在MATLAB命令窗中输入： syms x; y=sin(x)./x; eval(int(y,x,0,1) %eval()函数的功能就是将括号内的字符串视为语句并运行

11、三、结果：ans = 0.9461（3）;一、流程图：建立函数文件输入MATLAB指令 若不成立 修改结果二、源程序代码：建立被积函数文件funsfuns.mfunction f=funs(x)f=x.(-x);计算积分，在MATLAB命令窗中输入： quad(funs,0,1)三、结果：ans =1.2913（4）;一、流程图：建立函数文件输入MATLAB指令 若不成立 修改结果二、源程序代码：建立被积函数文件funffunf.mfunction f=funf(x)f=exp(2x)*sin(x).2;计算积分，在MATLAB命令窗中输入： d=pi/1000; x=0:d:2*pi; nt=length(x); y=funs