《Matlab 程序设计实践》课程考核.docx
《《Matlab 程序设计实践》课程考核.docx》由会员分享，可在线阅读，更多相关《《Matlab 程序设计实践》课程考核.docx（22页珍藏版）》请在冰豆网上搜索。
《Matlab程序设计实践》课程考核
《Matlab程序设计实践》课程考核
班级：
材料5562
学号：
054633898
姓名：
原迅
一．编程实现以下科学计算算法，并举一例应用之。
“帕德逼近”
在matlab中编程实现的帕德逼近法函数为：
Pade。
功能：
用帕德形式的有理分式逼近已知函数。
调用格式：
f=Pade（y,n）或f=Pade（y,n,x0）
其中，y为已知函数；
n为帕德有理分式的分母多项式的最高次数；
x0为逼近点的x坐标；
f为求得的帕德有理分式或在x0处的逼近值。
在matlab中实现函数的帕德逼近的代码如下：
M文件名：
Pade.m
functionf=pade（y,n,x0）
%用帕德形式的有理分式逼近已知函数
%已知函数：
y
%帕德有理分式的分母多项式的最高次数：
n
%逼近点的x坐标：
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=0;
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
方法一算法实现流程图如下：
No
Yes
No
Yes
No
Yes
No
Yes
No
Yes
No
Yes
No
Yes
No
Yes
例子：
帕德逼近应用实例，用勒让德公式（取四项）逼近函数1/（1-x），并求当x=0.5时的函数值
解：
在matlab命令窗口中输入以下命令：
M文件名：
aa.m
>>f=Pade（'1/（1-x）',4）
f=
（2.92857*t^4+0.821429*t^3+0.988095*t^2+1.0006*t+1.0）/（-0.5*t^4+0.107143*t^3-0.0119048*t^2+0.000595238*t+1.0）
M文件名：
bb.m
>>f=Pade（'1/（1-x）',4,0.5）
f=
2.0757
从逼近结果看，函数的准确值为1/（1-0.5）=2，而用4次有理分式的帕德逼近已经达到了很高的精度了。
实现流程图为：
二．编程解决以下科学计算和工程实际问题。
拉弯合成部件的截面设计。
这一设计计算将归结为解一个三次代数方程，过去要用试凑法反复运算，本例显示了用matlab求解的简洁。
钻床立柱如图所示。
设P=15kN，许用拉应力[σ]=35Mpa，钻头轴与立柱轴距离为0.4m，试求立柱直径。
钻床受力图
这个三次代数方程可用matlab多项式求根的roots函数求解
M文件名：
dd.m
源程序:
P=input（'p='）,l=input（'l='）%输入力和偏心距
asigma=input（'[Ò]='）%输入许用拉应力
a=[asigma*pi/32,0,-P/8,-P*1];%求三次代数方程的导数向量
r=roots（a）;%求代数方程的根
d=r（find（imag（r）==0））%只取实根
结果：
>>dd
p=15000
P=
15000
l=0.4
l=
0.4000
[σ]=35e6
asigma=
35000000
d=
0.1219
流程图如下：
三．求解下列边值问题：
（1）x”=（-2/t）x’+（2/t.^2）x+（10cos（ln（t）））/t.^2,其中x
（1）=1,x（3）=3,输出t=1.5,2,2.5时x的值，并作x的图。
解：
调用dsolve函数，求解微分方程的符号解法
M文件名：
ff.m
>>x=dsolve（'D2x=（-2/t）*D1x+（2/（t^2））*x-（10*cos（log（t）））/（t^2）','x
（1）=1','x（3）=3','t'）
x=
1/13*t*（28*tan（1/2*log（3））^2+9*tan（1/2*log（3））+1）/（1+tan（1/2*log（3））^2）-9/13/t^2*（6*tan（1/2*log（3））^2+tan（1/2*log（3））+3）/（1+tan（1/2*log（3））^2）-（2*tan（1/2*log（t））-3+3*tan（1/2*log（t））^2）/（1+tan（1/2*log（t））^2）
>>ezplot（x）
>>legend（'x-t的函数'）
>>t=[1.5,2,2.5]
t=
1.50002.00002.5000
>>x=subs（sym（x）,findsym（x）,t）
x=
2.47762.83362.9391
实现流程图如下：
（2）y”+（1/t）y’+（1-1/（4t.^2））y=sprt（t）cost,其中y
（1）=1,y（6）=-0.5,输出t=2,3,4,5时y的值，并作y的图。
源程序为：
M文件名：
gg.m
y=dsolve（'D2y+（1/t）*Dy+（1-1/（4*t^2））*y=sqrt（t）*cos（t）','y
（1）=1','y（6）=-0.5'）
t=[0:
0.1:
10]
y=subs（sym（y）,findsym（y）,t）
plot（t,y,'-r'）
legend（'y-t的函数'）
y=dsolve（'D2y+（1/t）*Dy+（1-1/（4*t^2））*y=sqrt（t）*cos（t）','y
（1）=1','y（6）=-0.5'）
t=[2,3,4,5]
y=subs（sym（y）,findsym（y）,t）
结果为：
>>gg
y=
（2^（1/2）*sin（t）*（（2^（1/2）*pi^（1/2）*t^2）/8-（2^（1/2）*pi^（1/2））/16+（2^（1/2）*pi^（1/2）*cos（2*t））/16+（2^（1/2）*pi^（1/2）*t*sin（2*t））/8））/（pi^（1/2）*t^（1/2））-（cos（t）*（sin（2*t）/4-（t*cos（2*t））/2））/（2*t^（1/2））-（6^（1/2）*sin（t）*（48*cos
（1）+16*6^（1/2）*cos（6）+142*6^（1/2）*cos
（1）*sin（6）-2*6^（1/2）*cos（6）*sin
（1）-4*6^（1/2）*cos（6）*sin
（1）*sin
（2）+24*6^（1/2）*cos
（1）*sin（6）*sin（12）-4*6^（1/2）*cos
（1）*cos
（2）*cos（6）+24*6^（1/2）*cos
（1）*cos（6）*cos（12）+2*6^（1/2）*cos
（1）*cos（6）*sin
（2）-2*6^（1/2）*cos
（2）*cos（6）*sin
（1）-2*6^（1/2）*cos
（1）*cos（6）*sin（12）+2*6^（1/2）*cos
（1）*cos（12）*sin（6）））/（96*t^（1/2）*（cos
（1）*sin（6）-cos（6）*sin
（1）））+（6^（1/2）*cos（t）*（48*sin
（1）+16*6^（1/2）*sin（6）+140*6^（1/2）*sin
（1）*sin（6）+2*6^（1/2）*cos
（1）*sin
（2）*sin（6）-2*6^（1/2）*cos
（2）*sin
（1）*sin（6）-2*6^（1/2）*cos（6）*sin
（1）*sin（12）+2*6^（1/2）*cos（12）*sin
（1）*sin（6）-4*6^（1/2）*sin
（1）*sin
（2）*sin（6）+24*6^（1/2）*sin
（1）*sin（6）*sin（12）-4*6^（1/2）*cos
（1）*cos
（2）*sin（6）+24*6^（1/2）*cos（6）*cos（12）*sin
（1）））/（96*t^（1/2）*（cos
（1）*sin（6）-cos（6）*sin
（1）））
t=
Columns1through9
00.10000.20000.30000.40000.50000.60000.70000.8000
Columns10through18
0.90001.00001.10001.20001.30001.40001.50001.60001.7000
Columns19through27
1.80001.90002.00002.10002.20002.30002.40002.50002.6000
Columns28through36
2.70002.80002.90003.00003.10003.20003.30003.40003.5000
Columns37through45
3.60003.70003.80003.90004.00004.10004.20004.30004.4000
Columns46through54
4.50004.60004.70004.80004.90005.00005.10005.20005.3000
Columns55through63
5.40005.50005.60005.70005.80005.90006.00006.10006.2000
Columns64through72
6.30006.40006.50006.60006.70006.80006.90007.00007.1000
Columns73throug
升级会员 