2019年實驗七matlab求解級數(shù)有關(guān)計算

實驗七matlab求解級數(shù)有關(guān)計算1.級數(shù)的基本概念常數(shù)項級數(shù)：稱用加號將數(shù)列"“的項連成的式子a+a+ah a+a+ah fah—123為(常數(shù)項)無窮級數(shù)，簡記為n=1。稱級數(shù)n=1前n項構(gòu)成的和S=a+a+a+?…+an1二工akk=1收斂，其和為S。limS=S收斂，其和為S。為級數(shù)的部分和。若n”n ，則稱級數(shù)n=1Taylor級數(shù)：設(shè)函數(shù)f(x)在包含x二a的區(qū)域內(nèi)具有各階導(dǎo)數(shù)，則稱幕級數(shù)￡f(n)(a)(x-a)n=f(a)+f'(a)(x-a)+f⑵(a)(x-a)2+???+f(n)(a)(x-a)n+…n! 2! n!n=0為函數(shù)f(X)在x=a的Taylor級數(shù)，當(dāng)a=0時稱為Maclaurin(麥克勞林)級數(shù)。2.級數(shù)的MATLAB命令MATLAB中主要用symsum,taylor求級數(shù)的和及進行Taylor展開。symsum(s,v,a,b)表達式s關(guān)于變量v從a到b求和taylor(f,a,n)將函數(shù)f在a點展為n-1階Taylor多項式可以用helpsymsum,helptaylor查閱有關(guān)這些命令的詳細信息例1先用taylor命令觀測函數(shù)y=sinX的Maclaurin展開式的前幾項，例如觀測前6項,相應(yīng)的MATLAB代碼為：>>clear;symsx;>>taylor(sin(x),0,1)>>taylor(sin(x),0,2)>>taylor(sin(x),0,3)>>taylor(sin(x),0,4)>>taylor(sin(x),0,5)>>taylor(sin(x),0,6)結(jié)果為:ans=0ans=xans=xans=x-1/6*xA3ans=x-1/6*xA3ans=xT/6*x3+l/120*x5然后在同一坐標(biāo)系里作出函數(shù)y二SinX和它的Taylor展開式的前幾項構(gòu)成的多項式函x3 x3x5TOC\o"1-5"\h\zy=x,y=x- ,y=x- + ，?…,數(shù) 3! 3! 5! 的圖形，觀測這些多項式函數(shù)的圖形向y二sinx的圖形的逼近的情況。例如，在區(qū)間［0，兀］上作函數(shù)y二sinx與多項式函數(shù)x3 x3x5y—x,y—x— ,y—x— + -3! 3! 5!圖形的MATLAB代碼為：>>x=0:0.01:pi;yl=sin(x);y2=x;y3=x-x43/6;y4=x-x43/6+x45/120;>>plot(x,yl,x,y2,':',x,y3,':',x,y4,':')類似地，根據(jù)函數(shù)的Taylor類似地，根據(jù)函數(shù)的Taylor級數(shù)結(jié)果如圖3.1，x2 x4 x6cosx—1— + — +?…,xe(—^o,+^o),2! 4! 6!ex=1+x+X2+X3+…,xe(-g,+x),2! 3!x2 x3 x4ln(1+x)—x— + — +?…,xe(—1,1],2 3 4a(a—1)(1+x)a—1+ax+ 2!x2+?…,xe(-1,1).作圖觀測其展開式的前幾項多項式逼近原函數(shù)的情況。例2利用冪級數(shù)計算指數(shù)函數(shù)。指數(shù)函數(shù)可展開為冪級數(shù)x2x3 xnex—1+x+ + +?…+—+?…，xe(—g,+8)2! 3! n!其通項為xM/prod(l:n),因此用下列循環(huán)相加就可計算出這個級數(shù)>>x=input('x=');n=input('n=');y=1;%輸入原始數(shù)據(jù)，初始化y>>fori=1:ny=y+x^i/prod(1:i);end,vpa(y,10),%將通項循環(huán)相加，得y執(zhí)行此程序，分別帶入x=1,2,4,-4這四個數(shù)，取n=10,y的結(jié)果如下2.7l828l80l，7.388994709，54.443l0406，.967l957672e-l而用vpa(exp(1),10),vpa(exp(2),10),vpa(exp(4),10),vpa(exp(-4),10)命令可得"，。‘創(chuàng)‘丘4的⑴位精確有效數(shù)字為2.718281828,7.389056099,54.59815003,.1831563889e-1對照可知，用級數(shù)法計算的有效數(shù)字分別為8，4，2，0位。由此可以看出，這個程序雖然原理上正確，但不好用。對不同的x,精度差別很大。其他存在的問題有：這個程序不能用于x的元素群運算；當(dāng)x為負數(shù)時，它成為交錯級數(shù)，收斂很慢；此n2程序要做2次乘法，n很大時，乘法次數(shù)太多，計算速度很低；對不同的x,要取不同的n才能達到精度要求，因此n不應(yīng)由用戶輸入，應(yīng)該由軟件按精度要求來選。正對上面的四個問題，可以采用下面四種方法改進：(1)允許數(shù)組輸入，改進輸出顯示x=input('x=');n=input('n=');y=ones(size(x));%輸入原始數(shù)據(jù)，初始化yfori=1:ny=y+x.Ai/prod(1:i);%循環(huán)相力口s1=sprintf('%13.0f',i);s2=sprintf('%15.8f',y);%將結(jié)果變?yōu)樽址甦isp([s1,s2]) %顯示end,執(zhí)行此程序，輸入x=[124-4]，n=10，結(jié)果為12.000000003.000000005.00000000-3.0000000022.500000005.0000000013.000000005.0000000032.666666676.3333333323.66666667-5.6666666742.708333337.0000000034.333333335.0000000052.716666677.2666666742.86666667-3.5333333362.718055567.3555555648.555555562.1555555672.718253977.3809523851.80634921-1.0952381082.718278777.3873015953.431746030.5301587392.718281537.3887125254.15414462-0.19223986102.718281807.3889947154.443104060.09671958可以利用exp(-x)=l/exp(x)來避免交錯級數(shù)的計算；為了減少乘法次數(shù)，設(shè)一個中間變量z,它的初始值為z=ones(size(x))，把循環(huán)體中的計算與句改為y=y+z;z=x.*z/i;這樣，求得的z就是z=x.T/i!，于是每個循環(huán)只需做一次乘法，計算整個級數(shù)只需n次乘法。按這種計算，y的初始值改為y=zeros(size(x))(4)為了按精度選擇循環(huán)次數(shù)，不該使用for循環(huán)，而用while語句，它可以設(shè)置循環(huán)的條件語句，通?？捎脃+z-y>tol,tol是規(guī)定的允許誤差.只要相鄰的兩次y值之差大于tol,循環(huán)就繼續(xù)進行，直到小于tol為止.xexp(x)二(expG-))k當(dāng)x較大時,exp(x)仍能很快收斂，還可以利用關(guān)系式 k，令x1=x/k.k通常取大于x而最接近x的2的幕，例如x=100，就取k=128，可以保證x1的絕對值小于1,這時級數(shù)收斂得很快??從練習(xí)中可以看出,n取10時(即級數(shù)取10項)就能保證7位有效數(shù)，而exp(X1)128可以化成X=(…((exp(X1))2)2…)2,即exp(x1)的7次自乘，總共用17次乘法就可完成exp(100)=(…((exp(1°0/128))2)2…)2的計算，這既保證了精度，又提高了速度.例3 編寫任意函數(shù)展開為各階泰勒級數(shù)的程序，并顯示其誤差曲線.對于任意函數(shù)y=f(x),其泰勒展開式為f(x)=f(a)+f'(a)(x-a)+-凹(x-a)2h f-凹(x-a)n+R(x).2! n! n其中Rn(x)為余項，也就是泰勒展開式的誤差.MATLAB語句為>>fxs=input('輸入y=f(x)的表達式','s'); %輸入原始條件,fxs是字符串>>K=input('輸入泰勒級數(shù)展開式的階K');>>a=input('展開的位置a=');>>b=input('展開的區(qū)間半寬度b=');>>x=linspace(a-b,a+b);%構(gòu)成自變量數(shù)組,確定其長度和步長>>lx=length(x);dx=2*b/(lx-1);>>y=eval(fxs);%求出y的準(zhǔn)確值>>subplot(1,2,1),plot(x,y,'.'),holdon%y的準(zhǔn)確值用點線繪出%求出a點的一階導(dǎo)數(shù)，注意求導(dǎo)后數(shù)組長度減少1>>Dy=diff(y)/dx;Dya(1)=Dy(round(lx-1)/2);>>yt(1,:)=y(round(lx/2))+Dya(1)*(x-a); %求y的一階泰勒展開,繪圖>>plot(x,yt(1,:))>>fork=2:K>>Dy=diff(y,k)/(dx^k);Dya(k)=Dy(round(lx-k)/2);%求a點k階導(dǎo)數(shù)>>yt(k,:)=yt(k-1,:)+Dya(k)/prod(1:k)*(x-a).^k; %求y的k階導(dǎo)

>>plot(x,yt(k,:));%繪圖>>e(k,:)=y-yt(k,:);%求出yt的誤差>>end>>title([fxs,'的各階泰勒級數(shù)曲線']),%注意如何組成標(biāo)注的字符串>>grid,holdoff,subplot(1,2,2)>>fork=1:Kplot(x,e(k,:)),holdon,end%繪制誤差曲線>>title([fxs,'的各階泰勒級數(shù)誤差曲線']),grid,holdoff執(zhí)行此程序，輸入fxs=cos(x),K=5,a=0.5,b=2,所得曲線見圖3.2(又變?yōu)檎`差曲線).讀者可以改變其坐標(biāo)系范圍以仔細觀測最關(guān)心的部分,也可輸入其他函數(shù)做驗算,注意輸入函數(shù)應(yīng)符合元素群運算規(guī)則.藝丄.例4 計算級數(shù)n=in2的值，可用symsum命令,相應(yīng)的MATLAB代碼為:>>clear;symsk;>>simple(symsum(1/kA2,1,Inf))%simple求解最簡形式,Inf為無窮大藝丄n藝丄n4n=1兀490結(jié)果為:ans=1/6*piA2類似地可驗證藝1兀6n6 945n=1蘭1可以猜想有n=1n22,2=1,2,其中m2是正整數(shù)，請驗證.藝1=e注：可用公式n=0n! 來計算e的近似值。如果要精確到小數(shù)點后15位，相應(yīng)的MATLAB代碼為:

>>digits(20); %設(shè)置今后數(shù)值計算以20位相對精度進行>>a=1.0;kk=1.0; %賦初值>>forn=1:20,kk=kk/n;,a=a+kk;,end>>vpa(a,17) %以17位相對精度給出a的值結(jié)果為e沁2.718281828459045.例5(調(diào)和級數(shù)例5(調(diào)和級數(shù))自然數(shù)的倒數(shù)組成的數(shù)列111231n 稱為調(diào)和數(shù)列，由調(diào)和數(shù)藝1 藝-列構(gòu)成的級數(shù)n=1n稱為調(diào)和級數(shù)，我們把它的前n項部分和k=1k記為H(n)。計算n=12，…,10時H(n)和Inn,ln(n+!)的值，并計算它們的差C(n)=H(n)—lnn，c(n)=H(n)—ln(n+1)，相應(yīng)的MATLAB代碼為:>>H(1)=1;C(1)=1;>>forn=2:100,H(n)=H(n-1)+1/n;,%for為循環(huán)語句>>c(n)=H(n)-log(n+1);,C(n)=H(n)-log(n);,end注意觀測C(n)單調(diào)遞減、C(n)單調(diào)遞增，二者相互接近的現(xiàn)象。計算n=10m(m=1,2,…,6)時C(n)和c(n)的值，注意觀測C(n)單調(diào)遞減、c(n)單調(diào)遞增，二者趨于同一極限的現(xiàn)象。并求出這個常數(shù)C。極限C=lim(1+ +—+?…+——lnn)=0.5771?…ntr 23n稱為歐拉(Euler)常數(shù)，顯然，可以證明它是一個無理數(shù)。c(n)=H(n)—ln(n+1)<C(n)=H(n)—lnn,C(n)—C(n)—c(n)=ln(1+丄)n趨于o,故C(n),c(n)趨于同一個極限C習(xí)題16-71.用taylor命令觀測函數(shù)y=f(x)的Maclaurin展開式的前幾項，然后在同一坐標(biāo)系里作出函數(shù)y=f(X)和它的Taylor展開式的前幾項構(gòu)成的多項式函數(shù)的圖形，觀測這些多項式函數(shù)的圖形向y=f(X)的圖形的逼近的情況(1)f(x)=arcsinx(2)(1)f(x)=arcsinx(2)f(x)=arctanx(3)f(x)=ex2ex(4)f(x