《列主元高斯消去法、LU三角分解法、龙贝格算法、最小二乘法的Matlab程序及运行结果》由会员分享,可在线阅读,更多相关《列主元高斯消去法、LU三角分解法、龙贝格算法、最小二乘法的Matlab程序及运行结果(7页珍藏版)》请在金锄头文库上搜索。
1、1列主元高斯消去法M文件functionx=gauss(a,b)n=length(a);x=zeros(n,1);a=ab;fork=1:n-1max=k;fori=k+1:nifa(i,k)a(max,k)max=i;endendtemp=a(k,k:n+1);a(k,k:n+1)=a(max,k:n+1);a(max,k:n+1)=temp;fori=k+1:na(i,k)=-a(i,k)/a(k,k);a(i,k+1:n+1)=a(i,k+1:n+1)+a(i,k)*a(k,k+1:n+1);endendx(n,1)=a(n,n+1)/a(n,n);fori=n-1:-1:1sum=0;
2、forj=i+1:nsum=sum+x(j,1)*a(i,j);endx(i,1)=(a(i,n+1)-sum)/a(i,i);endMatlab运行结果CommandWindow034J-11;212a.=0341-11212b-l;2;3b=123gaussans-(a,b)1.1667-0.33330.50002LU三角分解法M文件functiony=LU(A,B);n=length(A);A=AB;fork=1:n-1;fori=k:n;if(abs(A(i,k)=max(abs(A(k:n,k)P(k)=i;temp=A(k,:);A(k,:)=A(i,:);A(i,:)=temp;
3、endendforj=k+1:n;A(j,k)=A(j,k)/A(k,k);A(j,k+1:n+1)=A(j,k+1:n+1)-A(j,k)*A(k,k+1:n+1);endendP(n)=n;L(1,1)=1;L(2:n,1)=A(2:n,1);L(1,2:n)=0;U(1,1)=A(1,1);U(2:n,1)=0;U(1,2:n)=A(1,2:n);fori=2:n;L(i,1:i-1)=A(i,1:i-1);L(i,i)=1;L(i,i+1:n)=0;U(i,1:i-1)=0;U(i,i:n)=A(i,i:n);endx(n)=A(n,n+1)/U(n,n);fork=n-1:-1:1x
4、(k)=A(k,n+1);forp=n:-1:k+1;x(k)=x(k)-U(k,p)*x(p);endx(k)=x(k)/U(k,k);endxLUPEndMatlab运行结果CommandWindowa=034;1-11;21.2a二0341-11212b=:l;2;3b=123LU(a,b)x=L1667-|_33330.50QDL二LOOOO000L000000.5000-0.5000L00003龙贝格(Romberg)算法M文件functiont=romberg(f,a,b,e)t=zeros(15,4);t(1,1)=(b-a)/2*(f(a)+f(b);fork=2:4sum=0
5、;fori=10(k-2)sum=sum+f(a+(2*i-1)*(b-a)/2A(k-1);endt(k,1)=0.5*t(k-1,1)+(b-a)/2A(k-1)*sum;fori=2:kt(k,i)=(4A(i-1)*t(k,i-1)-t(k-1,i-1)/(4A(i-1)-1);endendfork=5:15sum=0;fori=1:2A(k-2)sum=sum+f(a+(2*i-1)*(b-a)/2A(k-1);endt(k,1)=0.5*t(k-1,1)+(b-a)/2A(k-1)*sum;fori=2:4t(k,i)=(4A(i-1)*t(k,i-1)-t(k-1,i-1)/(4
6、A(i-1)-1);endifk6ifabs(t(k,4)-t(k-1,4)=15disp(溢出);endMatlab运行结果CommandWindowf=inlineCsin(x)/xf,f)f二Iniinefunction:f(x)二sxn(x)/xrojnberg(fj10(-100),1,5*10(-?)Iff案0.94608ans二0,9207000(LB3BB0.9461000,94450.9461CL946100.94570.9461(L94610.94610.946D0.9461CL94610.94610.94fil0.9461(L94610.94610.94fil0.461H
7、94610.946100(00000Q000Q00000o0Q000000000004最小二乘法M文件functiona,max,det=zuixiaoerchengfa(x,y,r)n=length(x);c=ones(n,r+1);fori=2:r+1forj=1:nC(j,i)=X(j)F-1);endendA=C*C;b=C*y;a=inv(A)*b;det=0;maX=0;fori=1:nsum=a(1);forj=2:r+1sum=sum+a(j)*x(i)A(j-1);endCC=abs(y(i)-sum);ifCCmaxmax=CC;enddet=det+ccA2;enddet=sqrt(det);Matlab运行结果CommandWindowA=246SA二2468B=2112840B=2112840zuixiaoerchengfa(A3Bj2)ans=-8.76QD4.67500.1875zuixiaoerchengfa(A,B,1)arts-12.50006.5500