1、帕德逼近算法MATLAB程序设计实践课程作业一、用MATLAB编程实现“帕德逼近”的科学计算算法,及举例应用。1)帕德逼近算法说明如下:帕德逼近是一种有理分式逼近,逼近公式如下:大量实验表明,当L+M为常数时,取L=M,帕德逼近精确度最好,而且速度最快。此时,分子与分母中的系数可通过以下方式求解。首先,求解线形方程Aq=b,得到()的值,其中,然后,通过下式求出的值。注意,函数的帕德逼近不一定存在。在MATLAB中编程实现的帕德逼近法函数为:Pade。功能:用帕德形式的有理分式逼近已知函数。调用格式:f=Pade(y,n)或f=Pade(y,n,x0)。其中,y为已知函数; n为帕德有理分式的
2、分母多项式的最高次数; x0为逼近点的x坐标; f为求得的帕德有理分式或在x0处的逼近值。 2)程序源代码如下:在m文件中编写实现函数的Pade逼近的代码如下: function f=Pade(y,n,x0)%用帕德形式的有理分式逼近已知函数%已知函数:y%帕德有理分式的分母多项式的最高次数:n%逼近点的坐标:x0%求得的帕德有理分式或在x0处的逼近值:fsyms t;A=zeros(n,n);q=zeros(n,1);p=zeros(n+1,1);b=zeros(n,1);yy=0;a(1:2*n)=0.0;for(i=1:2*n) yy=diff(sym(y),findsym(sym(y)
3、,n); a(i)=subs(sym(yy),findsym(sym(yy),0.0)/factorial(i);end;for(i=1:n) for(j=1:n) A(i,j)=a(i+j-1); end; b(i,1)=-a(n+i);end; q=Ab;p(1)=subs(sym(y),findsym(sym(y),0.0);for(i=1:n) p(i+1)=a(n)+q(i)*subs(sym(y),findsym(sym(y),0.0); for(j=2:i-1) p(i+1)=p(i+1)+q(j)*a(i-j); endend f_1=0;f_2=1;for(i=1:n+1)
4、f_1=f_1+p(i)*(t(i-1);end for(i=1:n) f_2=f_2+q(i)*(ti);end if(nargin=3) f=f_1/f_2; f=subs(f,t,x0);else f=f_1/f_2; f=vpa(f,6);end3)算法实现流程图如下:4)用 勒让德 公式(取4项)逼近函数,并求当x=0.5时的函数值。 源代码及运算结果如下: f=Pade(1/(1-x),4) f = (1.+1.00060*t+.988095*t2+.821429*t3+2.92857*t4)/(1.+.595238e-3*t-.119048e-1*t2+.107143*t3-.5
5、00000*t4) f=Pade(1/(1-x),4,0.5)f =2.0757实现流程图为:二、求解工程问题,计算立柱的直径。已知:P=15kN, =35Mpa,l=0.4m1)建模:钻床受力图如下: 如图所示,钻床立柱受到拉力P和弯矩M=Pl的作用,立柱受到的拉应力之和为P与M的共同作用,其最大值(max)应小于许用拉应力,即:(max)=P/F+M/W (1)其中,对于圆柱体而言,F=d2/4为立柱截面积,W=d3/32为抗弯系数,将FW带入(1)式后,可得:=4P/d2+32M/d3 (2)化简上式,可得立柱直径d应满足:d3/32-Pd/8-Pl0 (3)于是,该工程问题可以看作是单
6、变量三次代数方程的求根问题!2)算法实现fzero函数: fzero函数用于求解单变量非线性方程根据(3)式的关系,将立柱直径d的系数项依次算出:=35*106MPa,P=15000N,l=0.4m,并带入到关系式中:在命令窗口输入: x=fzero(35*106*pi)/32*x3-(15000/8)*x-15000*0.4,0 1)x =0.1219流程图为:三、用MATLAB的符号解法,求解常微分方程:1: x=dsolve(D2x=(-2/t)*D1x+(2/(t2)*x-(10*cos(ln(t)/(t2),x(1)=1,x(3)=3) x = 1/13*t*(28*tan(1/2*
7、log(3)2+1+9*tan(1/2*log(3)/(1+tan(1/2*log(3)2)-9/13/t2*(6*tan(1/2*log(3)2+3+tan(1/2*log(3)/(1+tan(1/2*log(3)2)-(3*tan(1/2*log(t)2+2*tan(1/2*log(t)-3)/(1+tan(1/2*log(t)2) t=1.5,2,2.5t = 1.5000 2.0000 2.5000 x=subs(sym(x),findsym(x),t)x =2.4776 2.8336 2.9391plot(t,x,-r)流程图如下: 2 y=dsolve(D2y+(1/t)*Dy+(
8、1-1/(4*t2)*y=sqrt(t)*cos(t),y(1)=1,y(6)=-0.5) y = 1/4/t(1/2)*sin(t)*(-240*cos(1)*sin(1)4+160*cos(1)*sin(1)6-210*sin(1)+1330*sin(1)3-2240*sin(1)5+1120*sin(1)7+90*cos(1)*sin(1)2-4+72*sin(1)2-192*sin(1)4+128*sin(1)6-2*6(1/2)*cos(1)-5*cos(1)/(5-20*sin(1)2+16*sin(1)4)/sin(1)-1/2/t(1/2)*cos(t)*(-112*sin(1
9、)4+80*sin(1)6-105*sin(1)*cos(1)-560*cos(1)*sin(1)5+560*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(1/2)/(5-20*sin(1)2+16*sin(1)4)+1/4*(t*cos(t)+sin(t)*t2-sin(t)/t(1/2) t=2,3,4.5t = 2.0000 3.0000 4.5000 y=subs(sym(y),findsym(y),t)y =0.9544 -0.2187 -2.7915plot(t,y,-r)流程图如下:
copyright@ 2008-2022 冰豆网网站版权所有
经营许可证编号:鄂ICP备2022015515号-1