计算方法(本)大作业.docxVIP

《计算方法》大作业

1、根据拉格朗日n次插值多项式，编写拉格朗日插值多项式的

MATLAB通用程序(10分)

答：functionf=language(x,y,x0)

symst;

if(length(x)==length(y))

n=length(x);

else

disp(‘x和y的维数不相等！’);

return;

end

f=0.0;

for(i=1:n)

l=y(i);

for(j=1:i-1)

l=l*(t-x(j))/(x(i)-x(j));

end;

for(j=i+1:n)

l=l*(t-x(j))/(x(i)-x(j));

end;

f=f+1;

simplify(f);

if(i==n)

if(nargin==3)

f=subs(f,‘t’,x0);

else

f=collect(f);

f=vpa(f,6);

end

end

end

2、编写复化梯形公式的MATLAB通用程序(5分)

编程如下：

functions=mtrap(f,a,b,n);

symsx;

h=(b-a)/n;

x=linspace(a,b,n+1);

y=feval(f,x);

s=0.5*h*(y(1)+2*sum(y(2:n))+y(n+1));

3、编写龙贝格积分法的MATLAB通用程序(10分)

function[T]=mromb(f,a,b,eps)

ifnargin4,eps=1e-6;end

h=b-a;

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

n=1;J=0;err=1;

while(erreps)

J=J+1;h=h/2;S=0;

fori=1:n

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

S=S+feval(f,x);

end

R(J+1,1)=R(J,1)/2+h*S;

fork=1:J

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

end

err=abs(R(J+1,J+1)-R(J+1,J));

n=2*n;

end

R;

T=R(J+1,J+1);

4、编写列主元消元法的MATLAB通用程序(15分)

function[x]=mgauss2(A,b);

n=length(b);

fork=1:(n-1)

[ap,p]=max(abs(A(k:n,k)));

p=p+k-1;

ifpk

A([k,p],:)=A([p,k],:);

b([k,p],:)=b([p,k],:);

end

m=A(k+1:n,k)/A(k,k);

A(k+1:n,k+1:n)=A(k+1:n,k+1:n)-m*A(k,k+1:n);

b(k+1:n)=b(k+1:n)-m*b(k);

A(k+1:n,k)=zeros(n-k,1);

end

x=zeros(n,1);x(n)=b(n)/A(n,n);

fork=n-1:-1:1

x(k)=(b(k)-A(k,k+1:n)*x(k+1:n))/A(k,k);

end

5、编写直接三角分解法的MATLAB通用程序(15分)

function[x,L,R]=jc(A,b)

n=length(b);

R=zeros(n);

L=eye(n);

R(1,1:n)=A(1,1:n);

L(2:n,1)=A(2:n,1)/R(1,1);

fork=2:n

R(k,k:n)=A(k,k:n)-L(k,1:k-1)*R(1:k-1,k:n);

L(k+1:n,k)=(A(k+1:n,k)-L(k+1:n,1:k-1)*R(1:k-1,k))/R(k,k);

end

y=zeros(n,1);y(1)=b(1);

fork=2:n

y(k)=b(k)-L(k,1:k-1)*y(1:k-1);

end

x=zeros(n,1);x(n)=y(n)/R(n,n);

fork=n-1:-1:1

x(k)=(y(k)-R(k,k+1:n)*x(k+1:n))/R(k,k);

end

6、编写雅可比迭代法的MATLAB通用程序(10分)

function[x,i]=jacobi(A,b,x