计算方法上机题答案
- 格式：docx
- 大小：602.49 KB
- 文档页数：19
2.用下列方法求方程e^x+10x-2=0的近似根，要求误差不超过5*10的负4次方，并比较计算量
（1）二分法
（局部，大图不太看得清，故后面两小题都用局部截图）
（2）迭代法
（3）牛顿法
顺序消元法
#include
#include
#include
int main()
{ int N=4,i,j,p,q,k;
double m;
double a[4][5];
double x1,x2,x3,x4;
for (i=0;i
for (j=0;j
scanf("%lf",&a[i][j]);
for(p=0;p
{ for(k=p+1;k
{
m=a[k][p]/a[p][p];
for(q=p;q
a[k][q]=a[k][q]-m*a[p][q];
}
}
x4=a[3][4]/a[3][3];
x3=(a[2][4]-x4*a[2][3])/a[2][2];
x2=(a[1][4]-x4*a[1][3]-x3*a[1][2])/a[1][1];
x1=(a[0][4]-x4*a[0][3]-x3*a[0][2]-x2*a[0][1])/a[0][0];
printf("%f,%f,%f,%f",x1,x2,x3,x4);
scanf("%lf",&a[i][j]); （这一步只是为了看到运行的结果）
}
运行结果
列主元消元法
function[x,det,flag]=Gauss(A,b)
[n,m]=size(A);nb=length(b);
flag='OK';det=1;x=zeros(n,1);
for k=1:n-1 max1=0;
for i=k:n
if abs(A(i,k))>max1
max1=abs(A(i,k));r=i;
end
end if max1<1e-10
flag='failure';return;
end
if r>k
for j=k:n
z=A(k,j);A(k,j)=A(r,j);A(r,j)=z;
end
z=b(k);b(k)=b(r);b(r)=z;det=-det;
end
for i=k+1:n
m=A(i,k)/A(k,k);
for j=k+1:n
A(i,j)=A(i,j)-m*A(k,j);
end
b(i)=b(i)-m*b(k);
end
det=det*A(k,k);
end
det=det*A(n,n)
if abs(A(n,n))<1e-10
flag='failure';return;
end
x(n)=b(n)/A(n,n);
for k=n-1:-1：1
for j=k+1:n
b(k)=b(k)-A(k,j)*x(j);
end
x(k)=b(k)/A(k,k);
end
运行结果：
雅可比迭代法
function y=jacobi(a,b,x0)
D=diag(diag(a));
U=-triu(a,1);
L=-tril(a,-1);
B=D\(L+U);
f=D\b;
y=B*x0+f;n=1;
while norm(y-x0)>1e-4
x0=y;
y=B*x0+f;n=n+1;
end
y
n
高斯赛德尔迭代法
function y=seidel(a,b,x0)
D=diag(diag(a));
U=-triu(a,1);
L=-tril(a,-1);
G=(D-L)\U;
f=(D-L)\b;
y=G*x0+f;n=1;
while norm(y-x0)>10^(-4)
x0=y;
y=G*x0+f;n=n+1;
end
y
n
SOR迭代法
function y=sor(a,b,w,x0)
D=diag(diag(a));
U=-triu(a,1);
L=-tril(a,-1);
lw=(D-w*L)\((1-w)*D+w*U);
f=(D-w*L)\b*w;
y=lw*x0+f;n=1;
while norm(y-x0)>10^(-4)
x0=y;
y=lw*x0+f;n=n+1;
end
y
n
1.分段线性插值：
function y=fdxx(x0,y0,x)
p=length(y0);n=length(x0);m=length(x);
for i=1:m
z=x(i);
for j=1:n-1
if z
break;
end
end
y(i)= y0(j)*(z-x0(j+1))/(x0(j)-x0(j+1))+y0(j+1)*(z-x0(j))/(x0(j+1)-x0(j));
fprintf('y(%d)=%f\nx1=%.3fy1=%.3f\nx2=%.3fy2=%.3f\n\n',i,y(i),x0(j),y0(j),x0(j+1),y0(j+1));
end
end
结果0.39404 0.38007 0.35693
2.分段二次插值：
function y=fdec(x0,y0,x)
p=length(y0);n=length(x0);m=length(x);
for i=1:m
z=x(i);
for j=1:n-1
if z
break;
end
end
if j<2
j=j+1;
elseif (j
if (abs(x0(j)-z)>abs(x0(j+1)-z))
j=j+1;
elseif ((abs(x0(j)-z)==abs(x0(j+1)-z))&&(abs(x0(j-1)-z)>abs(x0(j+2)-z)))
j=j+1; end
end
ans=0.0;
for t=j-1:j+1
a=1.0;
for k=j-1:j+1
if t~=k
a=a*(z-x0(k))/(x0(t)-x0(k));
end
end
ans=ans+a*y0(t);
end
y(i)=ans;
fprintf('y(%d)=%f\n x1=%.3f y1=%.3f\n x2=%.3f y2=%.3f\n x3=%.3f
y3=%.3f\n\n',i,y(i),x0(j-1),y0(j-1),x0(j),y0(j),x0(j+1),y0(j+1));
end
end
结果为0.39447 0.38022 0.35725
3.拉格朗日全区间插值
function y=lglr(x0,y0,x)
p=length(y0);n=length(x0);m=length(x);
for t=1:m
ans=0.0;
z=x(t);
for k=1:n
p=1.0;
for q=1:n
if q~=k
p=p*(z-x0(q))/(x0(k)-x0(q));
end
end
ans=ans+p*y0(k);
end
y(t)=ans;
fprintf('y(%d)=%f\n',t,y(t));
end
end
结果为0.39447 0.38022 0.35722
