2、flag='failure'return;endfork=i+1:nm=0;forq=1:i-1m=m+L(k,q)*U(q,i);endL(k,i)=(A(k,i)-m)/U(i,i);endendLUend1.結(jié)果>>LU_decom(A)m=4n=4L=1000210012103321U=24260123003600015.拓展在編寫程序過程中由于角標(biāo)較多因此在運(yùn)行過程中出現(xiàn)了不少角標(biāo)不對的錯誤題目;給出函數(shù)f(x)=1/(1+25x^2),求f(x在[-1,1]上取5個和9個等距節(jié)點(diǎn),做最小二乘擬
3、合,得出均方誤差方誤差。五個節(jié)點(diǎn)時,matlab編碼為:首先建立M文件,并保存functiony=f(x)y=1/(1+25*x^2);endx=[-1-0.500.51];fori=1:5y(i)=f(x(i));enda=polyfit(x,y,3)symsxf1=a(1)*x^3+a(2)*x^2+a(3)*x+a(4)x=[-1-0.500.51];fori=1:5E(i)=(f(x(i))-(a(1)*x(i)^3+a(2)*x(i)^2+a(3)*x(i)+a(4)))^2;endsum(E)
4、輸出結(jié)果為a=-0.0000-0.6063-0.00000.5737f1=-47781/4*x^3-1600/2639*x^2-30703/*x+1514/2639(擬合的多項式)ans=0.3534(均方誤差)九個點(diǎn)的時候,matlab編碼為:x=[-1-0.75-0.5-0.2500.250.50.751];fori=1:9y(i)=f(x(i));enda=polyfit(x,y,3)symsxf2=a(1)*x^3+a(2)*x^2+a(3)*x+a(4)x=[-1-0.75-0.5-0.2500
5、.250.50.751];fori=1:5E1(i)=(f(x(i))-(a(1)*x(i)^3+a(2)*x(i)^2+a(3)*x(i)+a(4)))^2;endsum(E1)輸出結(jié)果為:a=-0.0000-0.56090.00000.4855f2=-6039/*x^3-12921/85248*x^2+54361/*x+64671/85248(最小二乘擬合多項式)ans=0.3350(均方誤差)用復(fù)合梯形公式求積分的值。functioni=combinetraprl(f,a,b,eps)%復(fù)化梯形公式
6、求函數(shù)f在區(qū)間[a,b]上的定積分%函數(shù)名:f%積分下限:a%積分上限:b%積分精度:eps%積分值:i%積分劃分的子區(qū)間個數(shù):stepn=1;h=(b-a)/2;i1=0;i2=(subs(sym(f),findsym(sym(f)),a)+subs(sym(f),findsym(sym(f)),b))/h;whileabs(i2-i1)>epsn=n+1;h=(b-a)/n;i1=i2;i2=0;fori=0:n-1x=a+h*i;x1=x+h;i2=i2+(h/2)*(subs(sym(f),fin
7、dsym(sym(f)),x)+...subs(sym(f),findsym(sym(f)),x1));endend四階龍格-庫塔法分別求解下列初值問題;function[x,y]=runge_kutta(fun,x0,xt,y0,pointnum,varargin)ifnargin<3y0=0;endy(1,:)=y0(:)';h=(xt-x0)/(pointnum-1);x=x0+[0:pointnum]'*h;fork=1:pointnumf1=h*feval(fun,x(k),y(k,:),var
8、argin{:});f1=f1(:)';f2=h*feval(fun,x(k)+h/2,y(k,:)+f1/2,varargin{:});f2=f2(:)';f3=h*feval(fun,x(k)+h/2,y(k,:)+f2/2,varargin{:});f3=f3(:)';f4=h*feval(fun,x(k)+h,y(k,:)+f3,varargin{:});f4=f4(:)';y(k+1,:)=y(k,:)+(f1+2*