帕德逼近算法.docx
- 文档编号:30535338
- 上传时间:2023-08-16
- 格式:DOCX
- 页数:15
- 大小:162.25KB
帕德逼近算法.docx
《帕德逼近算法.docx》由会员分享,可在线阅读,更多相关《帕德逼近算法.docx(15页珍藏版)》请在冰豆网上搜索。
帕德逼近算法
《MATLAB程序设计实践》课程作业
一、用MATLAB编程实现“帕德逼近”的科学计算算法,及举例应用。
1)帕德逼近算法说明如下:
帕德逼近是一种有理分式逼近,逼近公式如下:
大量实验表明,当L+M为常数时,取L=M,帕德逼近精确度最好,而且速度最快。
此时,分子与分母中的系数可通过以下方式求解。
首先,求解线形方程Aq=b,得到(
…
)的值,其中
,
,
然后,通过下式求出
的值。
注意,函数的帕德逼近不一定存在。
在MATLAB中编程实现的帕德逼近法函数为:
Pade。
功能:
用帕德形式的有理分式逼近已知函数。
调用格式:
f=Pade(y,n)或f=Pade(y,n,x0)。
其中,y为已知函数;
n为帕德有理分式的分母多项式的最高次数;
x0为逼近点的x坐标;
f为求得的帕德有理分式或在x0处的逼近值。
2)程序源代码如下:
①在m文件中编写实现函数的Pade逼近的代码如下:
functionf=Pade(y,n,x0)
%用帕德形式的有理分式逼近已知函数
%已知函数:
y
%帕德有理分式的分母多项式的最高次数:
n
%逼近点的坐标:
x0
%求得的帕德有理分式或在x0处的逼近值:
f
symst;
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)),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=A\b;
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);
end
end
f_1=0;
f_2=1;
for(i=1:
n+1)
f_1=f_1+p(i)*(t^(i-1));
end
for(i=1:
n)
f_2=f_2+q(i)*(t^i);
end
if(nargin==3)
f=f_1/f_2;
f=subs(f,'t',x0);
else
f=f_1/f_2;
f=vpa(f,6);
end
3)算法实现流程图如下:
4)用勒让德公式(取4项)逼近函数
,并求当x=0.5时的函数值。
源代码及运算结果如下:
>>f=Pade('1/(1-x)',4)
f=
(1.+1.00060*t+.988095*t^2+.821429*t^3+2.92857*t^4)/(1.+.595238e-3*t-.119048e-1*t^2+.107143*t^3-.500000*t^4)
>>f=Pade('1/(1-x)',4,0.5)
f=
2.0757
实现流程图为:
二、求解工程问题,计算立柱的直径。
已知:
P=15kN,[σ]=35Mpa,l=0.4m
1)建模:
钻床受力图如下:
如图所示,钻床立柱受到拉力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-Pl≥0(3)
于是,该工程问题可以看作是单变量三次代数方程的求根问题!
2)算法实现
fzero函数:
fzero函数用于求解单变量非线性方程
根据(3)式的关系,将立柱直径d的系数项依次算出:
[σ]=35*106MPa,P=15000N,l=0.4m,并带入到关系式中:
在命令窗口输入:
>>x=fzero('(35*10^6*pi)/32*x^3-(15000/8)*x-15000*0.4',[01])
x=
0.1219
流程图为:
三、用MATLAB的符号解法,求解常微分方程:
1:
>>x=dsolve('D2x=(-2/t)*D1x+(2/(t^2))*x-(10*cos(ln(t)))/(t^2)','x
(1)=1','x(3)=3')
x=
1/13*t*(28*tan(1/2*log(3))^2+1+9*tan(1/2*log(3)))/(1+tan(1/2*log(3))^2)-9/13/t^2*(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.5]
t=
1.50002.00002.5000
>>x=subs(sym(x),findsym(x),t)
x=
2.47762.83362.9391
plot(t,x,'-r')
流程图如下:
2>>y=dsolve('D2y+(1/t)*Dy+(1-1/(4*t^2))*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)^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)*t^2-sin(t))/t^(1/2)
>>t=[2,3,4.5]
t=
2.00003.00004.5000
>>y=subs(sym(y),findsym(y),t)
y=
0.9544-0.2187-2.7915
plot(t,y,'-r')
流程图如下:
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 逼近 算法
![提示](https://static.bdocx.com/images/bang_tan.gif)