帕德逼近算法.docx
MATLAB程序设计实践课程作业一、用MATLAB编程实现“帕德逼近的科学计算算法，及举例应用。13帕德逼近算法说明如下：帕德逼近是一种有理分式逼近，逼近公式如下：大量实验说明，当L+M为常数时，取L=M,帕德逼近精确度最好，而且速度最快。此时，分子与分母中的系数可通过以下方式求解。首先，求解线形方程Aq=b,得到(%,%,%/)的值，其中«1a2an或4+1-lJ|_一。2,1然后，通过下式求出(Po,P,P2P”)的值。注意，函数的帕德逼近不一定存在。在MATLAB中编程实现的帕德逼近法函数为：Pade0功能：用帕德形式的有理分式逼近函数。调用格式:f=Pacle(y,n)或f=Pade(y,n,x)。其中，y为函数；n为帕德有理分式的分母多项式的最高次数;为逼近点的X坐标；f为求得的帕德有理分式或在x处的逼近值。2程序源代码如下：在m文件中编写实现函数的Pade逼近的代码如下:functionf=Pade(y,n,x)%用帕德形式的有理分式逼近函数%函数：y%帕德有理分式的分母多项式的最高次数：n%逼近点的坐标：x%求得的帕德有理分式或在x处的逼近值：fsymst;A=zeros(n,n);q=zeros(n,1);p=zeros(n+l,1);b=zeros(n,1);yy=O;a(l:2*n)=0.O;for(i=l:2*n)yy=diff(sym(y),findsym(sym(y),n);a(i)=subs(sym(yy),findsym(sym(yy),O.O)/factorial(i);end;for(i=l:n)for(j=l:n)A(i,j)=a(i+j-l);end;b(i,l)=-a(n+i);end;q=Ab;p(l)=subs(sym(y),findsym(sym(y),O.O);for(i=l:n)p(i+l)=a(n)+q(i)*subs(sym(y),findsym(sym(y),O.O);for(j=2:i-l)p(i+l)=p(i+l)+q(j)*a(i-j);endendf_l=O;f_2=l;for(i=kn+l)f-l=f-l÷p(i)*(t(i-D);endfor(i=l:n)f_2=f_2+q(i)*(ti);endif(nargin-3)f=f_l/f_2;f=subs(f,t,x);elsef=f_l/f_2;f=vpa(f,6);end3算法实现流程图如下：开始/定义变量，输入:SymS t;/A=zeros(n,n);q=zeros(n, 1);/p=zeros(n+1,1 );b=zeros(n, 1);源代码及运算结果如下:>>f=Pade('l(l-x),4)f=(1.+1.00060*t+.988095*t2+.821429*t3+2.92857*t4)(1.+.595238e-3*t-.119048e-l*t2+.107143*tA3-.5OOOOO*tA4)>>f=Pade(,l(l-x)',4,0.5)f=2.0757实现流程图为:二、求解工程问题，计算立柱的直径。:P=15kN,=35Mpa,l=0.4mU建模：钻床受力图如下：I如下图，钻床立柱受到拉力P和弯矩M=PI的作用，立柱受到的拉应力之和为P与M的共同作用，其最大值。(max)应小于许用拉应力。,即：(ma)=PF+MW(1)其中，对于圆柱体fL一一工m,_抗弯系数，将FW带入(1)式后，可得：=4Pd2÷32Md3。J(2)化简上d332-Pd8-Pl0(3)于是，2算法实现fzero函数：fzero函数用于求解单变量非线性方程根据式的关系，将立柱直径d的系数项依次算出：。=35*106MPa,P=15000N,l=0.4m,并带入到关系式中：在命令窗口输入：»x=fzero('(35*106*pi)32*x3-(150008)*x-15000*0.401)X=0.1219流程图为：三、用MATLAB的符号解法,求解常微分方程：1：»x=dsolve('D2x=(-2t)*D1x+(2(t2)*x-(10*cos(ln(t)(t2),'x(1)=1'x(3)=3')X=113*t*(28*tan(l2*log(3)2+1+9*tan(12*log(3)(1+tan(12*log(3)2)-913t2*(6*tan(l2*log(3)2+3+tan(12*log(3)(l+tan(l2*log(3)2)-(3*tan(l2*log(t)2+2*tan(l2*log(t)-3)(l+tan(l2*log(t)2)»t=l,5,2,2.5t=1.50002.00002.5000»X=SUbS(Sym(X),findsym(x),DX=2.47762.83362.9391PIOt(Lx,'-r')流程图如下：cos(1)*sin(1)3+35*sin(1)2-12*cos(1)-64*cos(1)*sin(1)4+64*cos(1)*sin(1)2-6(12)(5-20*sin(1)2+16*sin(1)4)+14*(t*cos(t)+sin(t)*t2-sin(t)t(l/2)»t=2,3,4.5Jt=2.00003.00004.5000»y=subs(sym(y),findsym(y),t)y=0.9544-0.2187-2.7915plot(t,y,'-r')流程图如下：
