计算方法中涉及的各种算法程序与结果

列主元消去法function [X]=Gauss(A,B)

[m,n]=size(A);

l=length(B);

if n~=m

fprintf('wrong');

return;end

if m~=l

fprintf('wrong');

return;end

fprintf('go to calculate');

for k=1:n-1

max=0;

for j=k:n

if abs(A(j,k))>max

max=abs(A(j,k));mark=j;

end

end

if abs(max)<1e-5

fprintf('wrong');

return;end

for i=k:n

a=A(k,i);A(k,i)=A(mark,i);A(mark,i)=a;

end

b=B(k);B(k)=B(mark);B(mark)=b;

for j=k+1:n

m=A(j,k)/A(k,k);

for i=k:n

A(j,i)=A(j,i)-m*A(k,i);

end

B(j)=B(j)-m*B(k);

end

End

for j=n:-1:1

for i=j+1:n

B(j)=B(j)-A(j,i)*X(i);

end

X(j)=B(j)/A(j,j);

end

三角分解法function [L,U,y,x]=LU(A,B)

n=length(A);

L=eye(n);

U=zeros(n,n);

for k=1:n

for j=k:n

c=0;

for i=1:k-1

c=c+L(k,i)*U(i,j);

end

U(k,j)=A(k,j)-c;

end

if U(k,k)==0

for i=k+1:n

z=i;

if U(k,i)~=0

break;end

end

for i=1:k

a=U(i,k);U(i,k)=U(i,z);U(i,z)=a;

end

for i=1:n

a=A(i,k);A(i,k)=A(i,z);A(i,z)=a;

end

end

for i=k+1:n

c=0;

for q=1:k-1

c=c+L(i,q)*U(q,k);

end

L(i,k)=(A(i,k)-c)/U(k,k);

end

end

for k=1:n

c=0;

for i=1:k-1

c=c+L(k,i)*y(i);

end

y(k)=(B(k)-c)/L(k,k);

end

for k=n:-1:1

c=0;

for i=k+1:n

c=c+U(k,i)*x(i);

end

x(k)=(y(k)-c)/U(k,k); end

最小二乘法拟合多项式function y = curve(x,y,n)

L=size(x);

m=L(2);

for i=1:m;

D(i,1)=1;

for j=2:n+1;

D(i,j)=D(i,j-1)*x(i);

end

end

A=D'*D;

b=D'*y';

S=A\b;

for i=1:n+1

R(i)=S(n+2-i);

end

y=poly2str(R,'x');

YY=D*S;

sum=0;

max=0;

for i=1:m

D(i)=abs(YY(i)-y(i));

if max

max=D(i);

end

sum=sum+D(i)^2;

end

sum=sqrt(sum);

disp(y);

disp('均方误差为');

disp(sum);

disp('最大误差为');

disp(max);

end

龙贝格

function [Y,R]=rom(a,b,e,f)

m=1;

d=1;

h=b-a;

i=0;

R=zeros(4,4);

R(1,1)=h*(f(a)+f(b))/2;

while d>e

i=i+1;

h=h/2;

s=0;

for p=1:m

x=a+h*(2*p-1);

s=s+f(x);

end

R(i+1,1)=R(i,1)/2+h*s;

m=m*2;

p=i;

if p>3

p=3;

end

for k=1:p

R(i+1,k+1)=((4^k)*R(i+1,k)-R(i,k))/(4^k-1);

end

d=abs(R(i,p)-R(i+1,p));

end

Y=R(i+1,4);