《MATLAB程序設(shè)計(jì)實(shí)踐》課程作業(yè)一、用MATLAB編程實(shí)現(xiàn)“帕德逼近〞的科學(xué)計(jì)算算法，及舉例應(yīng)用。1〕帕德逼近算法說(shuō)明如下：帕德逼近是一種有理分式逼近，逼近公式如下：大量實(shí)驗(yàn)說(shuō)明，當(dāng)L+M為常數(shù)時(shí)，取L=M，帕德逼近精確度最好，而且速度最快。此時(shí)，分子與分母中的系數(shù)可通過(guò)以下方式求解。首先，求解線形方程Aq=b，得到〔…〕的值，其中，，然后，通過(guò)下式求出的值。注意，函數(shù)的帕德逼近不一定存在。在MATLAB中編程實(shí)現(xiàn)的帕德逼近法函數(shù)為：Pade。功能：用帕德形式的有理分式逼近函數(shù)。調(diào)用格式：f=Pade(y,n)或f=Pade(y,n,x0)。其中，y為函數(shù)；n為帕德有理分式的分母多項(xiàng)式的最高次數(shù)；x0為逼近點(diǎn)的x坐標(biāo)；f為求得的帕德有理分式或在x0處的逼近值。2〕程序源代碼如下：①在m文件中編寫(xiě)實(shí)現(xiàn)函數(shù)的Pade逼近的代碼如下：functionf=Pade(y,n,x0)%用帕德形式的有理分式逼近函數(shù)%函數(shù)：y%帕德有理分式的分母多項(xiàng)式的最高次數(shù)：n%逼近點(diǎn)的坐標(biāo)：x0%求得的帕德有理分式或在x0處的逼近值：fsymst;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);endendf_1=0;f_2=1;for(i=1:n+1)f_1=f_1+p(i)*(t^(i-1));endfor(i=1:n)f_2=f_2+q(i)*(t^i);endif(nargin==3)f=f_1/f_2;f=subs(f,'t',x0);elsef=f_1/f_2;f=vpa(f,6);end3〕算法實(shí)現(xiàn)流程圖如下：開(kāi)始開(kāi)始定義變量，輸入：定義變量，輸入：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.0No開(kāi)始循環(huán)判斷No開(kāi)始循環(huán)判斷i≤2nYesYesyy=diff(sym(y),findsym(sym(y)),n);yy=diff(sym(y),findsym(sym(y)),n);a(i)=subs(sym(yy),findsym(sym(yy)),0.0)/factorial(i);開(kāi)始循環(huán)判斷開(kāi)始循環(huán)判斷i≤nNoNoYesYesNo開(kāi)始循環(huán)判斷No開(kāi)始循環(huán)判斷j≤nq=A\b;q=A\b;p(1)=subs(sym(y),findsym(sym(y)),0.0);b(i,1)=-a(n+i)Yesb(i,1)=-a(n+i)YesA(i,j)=a(i+j-1)A(i,j)=a(i+j-1)No開(kāi)始循環(huán)判斷No開(kāi)始循環(huán)判斷j≤nYesYesp(i+1)=a(n)+q(i)*subs(sym(y),findsym(sym(y)),0.0)p(i+1)=a(n)+q(i)*subs(sym(y),findsym(sym(y)),0.0)f_1=0;f_1=0;f_2=1;No開(kāi)始循環(huán)判斷No開(kāi)始循環(huán)判斷2≤j≤i-1YesYesp(i+1)=p(i+1)+q(j)*a(i-j)p(i+1)=p(i+1)+q(j)*a(i-j)開(kāi)始循環(huán)判斷開(kāi)始循環(huán)判斷1≤j≤n+1NoNoYesYesf_1=f_1+p(i)*(t^(i-1))f_1=f_1+p(i)*(t^(i-1))開(kāi)始循環(huán)判斷1開(kāi)始循環(huán)判斷1≤i≤nNoNoYesYesf_2=f_2+q(i)*(t^i)f_2=f_2+q(i)*(t^i)對(duì)對(duì)nargin==3進(jìn)行循環(huán)判斷NoNoYesYes輸出結(jié)果f=f_1/f_2;輸出結(jié)果f=f_1/f_2;f=vpa(f,6)輸出結(jié)果f=f_1/f_2;f=subs(f,'t',x0);結(jié)束結(jié)束4〕用勒讓德公式〔取4項(xiàng)〕逼近函數(shù)，并求當(dāng)x=0.5時(shí)的函數(shù)值。源代碼及運(yùn)算結(jié)果如下：>>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實(shí)現(xiàn)流程圖為：開(kāi)始開(kāi)始輸入：f=Pade('1/(1-x)',4)輸出結(jié)果，其形式為帕德有理分式調(diào)用帕德函數(shù)，分母最高次項(xiàng)為4次調(diào)用帕德函數(shù)，分母最高次項(xiàng)為4次，計(jì)算x=0.5時(shí)的逼近值輸入：f=Pade(‘1/(1-x)’,4,0.5)輸出結(jié)果，其形式為x=0.5時(shí)的逼近值結(jié)束二、求解工程問(wèn)題，計(jì)算立柱的直徑。：P=15kN,[σ]=35Mpa，l=0.4m1〕建模：鉆床受力圖如下：llPPM=PlA鉆床受力圖AA立柱受力圖A-A如下圖，鉆床立柱受到拉力P和彎矩M=Pl的作用，立柱受到的拉應(yīng)力之和為P與M的共同作用，其最大值σ〔max〕應(yīng)小于許用拉應(yīng)力[σ]，即：σσ(max)=P/F+M/W≤[σ]〔1〕其中，對(duì)于圓柱體而言，F(xiàn)=πd2/4為立柱截面積，W=πd3/32為抗彎系數(shù)，將FW帶入〔1〕式后，可得：σσ=4P/πd2+32M/πd3≤[σ](2)化簡(jiǎn)上式，可得立柱直徑d應(yīng)滿足：[[σ]πd3/32-Pd/8-Pl≥0〔3〕于是，該工程問(wèn)題可以看作是單變量三次代數(shù)方程的求根問(wèn)題！2〕算法實(shí)現(xiàn)fzero函數(shù):fzero函數(shù)用于求解單變量非線性方程根據(jù)〔3〕式的關(guān)系，將立柱直徑d的系數(shù)項(xiàng)依次算出：[σ]=35*106MPa，P=15000N，l=0.4m，并帶入到關(guān)系式中：在命令窗口輸入：>>x=fzero('(35*10^6*pi)/32*x^3-(15000/8)*x-15000*0.4',[01])x=0.1219流程圖為：開(kāi)始開(kāi)始調(diào)用fzero函數(shù)輸入：x=fzero(f，x0)輸出：區(qū)間[0，1]的根，x=0.1219m結(jié)束三、用MATLAB的符號(hào)解法,求解常微分方程: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.9391plot(t,x,'-r')流程圖如下：x=x=dsolve('D2x=(-2/t)*D1x+(2/(t^2))*x-(10*cos(ln(t)))/(t^2)''x(1)=1','x(3)=3''x(1)=1','x(3)=3'xx=?t=[1.5,2,2.5]x=subs(sym(x),findsym(x),t)t=[1.5,2,2.5]x=subs(sym(x),findsym(x),t)x=x=?plot(t,x,'-r')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.7915plot(t,y,'-r')流程圖如下：y=dsolve('D2y+(1/t)*Dy+(1-1/(4*t^2))*y=sqrt(t)