《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

function f= pade(y,n,x0)

%用帕德形式的有理分式逼近已知函数

%已知函数：y

%帕德有理分式的分母多项式的最高次数：n

%逼近点的x坐标：x0

%求得的帕德有理分式或在x0处的逼近值：f

syms 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)),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

方法一算法实现流程图如下：

例子：帕德逼近应用实例，用勒让德公式（取四项）逼近函数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.5000

2.0000 2.5000

>> x=subs(sym(x),findsym(x),t)

x =

2.4776 2.8336 2.9391

1/13 t (28 tan(1/2 log(3))2+9 tan(1/2 log(3))+1)/(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)

t