帕德逼近算法.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'）
流程图如下：
升级会员 