计算方法上机作业集合1.docx
《计算方法上机作业集合1.docx》由会员分享,可在线阅读,更多相关《计算方法上机作业集合1.docx(31页珍藏版)》请在冰豆网上搜索。
计算方法上机作业集合1
第一次&第二次上机作业
上机作业:
1.在Matlab上执行:
>>5.1-5-0.1和>>1.5-1-0.5
给出执行结果,并简要分析一下产生现象的原因。
解:
执行结果如下:
在Matlab中,小数值很难用二进制进行描述。
由于计算精度的影响,相近两数相减会出现误差。
2.(课本181页第一题)
解:
(1)n=0时,积分得
=ln6-ln5,编写如下图代码
从以上代码显示的结果可以看出,
的近似值为0.012712966517465
(2)
=
dx,可得
dx
dx
dx,得
则有
取
以此逆序估算
。
代码段及结果如下图所示
结果是从
逆序输出至
,所以得到
的近似值为.0883********。
(3)从
估计的过程更为可靠。
首先根据积分得表达式是可知,被积函数随着n的增大,其所围面积应当是逐步减小的,即积分值应是随着n的递增二单调减小的,
(1)中输出的值不满足这一条件,
(2)满足。
设
表示
的近似值,
-
=
(根据递推公式可以导出此式),可以看出,随着n的增大,误差也在增大,所以顺序估计时,算法不稳定性逐渐增大,逆序估计情况则刚好相反,误差不断减小,算法逐渐趋于稳定。
2.(课本181页第二题)
(1)上机代码如图所示
求得近似根为0.09058
(2)上机代码如图所示
得近似根为0.09064;
(3)牛顿法上机代码如下
计算所得近似解为0.09091
第三次上机作业
上机作业181页第四题
线性方程组为
=
(1)顺序消元法
A=[1.1348,3.8326,1.1651,3.4017;0.5301,1.7875,2.5330,1.5435;3.4129,4.9317,8.7643,1.3142;1.2371,4.9998,10.6721,0.0147];
b=[9.5342;6.3941;18.4231;16.9237];
上机代码(函数部分)如下
function[b]=gaus(A,b)%用b返回方程组的解
B=[A,b];
n=length(b);
RA=rank(A);
RB=rank(B);
dif=RB-RA;
ifdif>0
disp('此方程组无解');
return
end
ifRA==RB
ifRA==n
formatlong;
disp('此方程组有唯一解');
forp=1:
n-1
fork=p+1:
n
m=B(k,p)/B(p,p);
B(k,p:
n+1)=B(k,p:
n+1)-m*B(p,p:
n+1);
end
end%顺序消元形成上三角矩阵
b=B(1:
n,n+1);
A=B(1:
n,1:
n);
b(n)=b(n)/A(n,n);
forq=n-1:
-1:
1
b(q)=(b(q)-sum(A(q,q+1:
n)*b(q+1:
n)))/A(q,q);
end%回代求解
else
disp('此方程组有无数组解');
end
end
上机运行结果为
>>A=[1.1348,3.8326,1.1651,3.4017;0.5301,1.7875,2.5330,1.5435;3.4129,4.9317,8.7643,1.3142;1.2371,4.9998,10.6721,0.0147];
b=[9.5342;6.3941;18.4231;16.9237];
>>X=gaus(A,b)
此方程组有唯一解
X=
1.000000000000000
1.000000000000000
1.000000000000000
1.000000000000000
>>
(2)列主元消元法
A=[1.1348,3.8326,1.1651,3.4017;0.5301,1.7875,2.5330,1.5435;3.4129,4.9317,8.7643,1.3142;1.2371,4.9998,10.6721,0.0147];
b=[9.5342;6.3941;18.4231;16.9237];
函数部分代码如下
function[b]=lieZhu(A,b)%用b返回方程组的解
B=[A,b];
RA=rank(A);
RB=rank(B);
n=length(b);
dif=RB-RA;
formatlong;
ifdif>0
disp('该方程组无解');
return
end
ifdif==0
ifRA==n
disp('该方程组有唯一解');
c=zeros(1,n);
fori=1:
n-1
max=abs(A(i,i));
m=i;
forj=i+1:
n
ifmaxmax=abs(A(j,i));
m=j;
end
end%求出每一次消元时绝对值最大的一行的行号
ifm~=i
fork=i:
n
c(k)=A(i,k);
A(i,k)=A(m,k);
A(m,k)=c(k);
end
d1=b(i);
b(i)=b(m);
b(m)=d1;%函数值跟随方程一起换位置
end
fork=i+1:
n
forj=i+1:
n
A(k,j)=A(k,j)-A(i,j)*A(k,i)/A(i,i);
end
b(k)=b(k)-b(i)*A(k,i)/A(i,i);
A(k,i)=0;
end
end%完成消元操作,形成上三角矩阵
b(n)=b(n)/A(n,n);
fori=n-1:
-1:
1
sum=0;
forj=i+1:
n
sum=sum+A(i,j)*b(j);
end
b(i)=(b(i)-sum)/A(i,i);%回代求解其他未知数
end
end
else
disp('此方程组有无数组解');
end
end上机运行,结果为
>>A=[1.1348,3.8326,1.1651,3.4017;0.5301,1.7875,2.5330,1.5435;3.4129,4.9317,8.7643,1.3142;1.2371,4.9998,10.6721,0.0147];
b=[9.5342;6.3941;18.4231;16.9237];
X=lieZhu(A,b)
该方程组有唯一解
X=
1.000000000000000
1.000000000000002
0.999999999999999
0.999999999999999
>>
根据两种方法运算结果,顺序消元法得到结果比列主元消元法得到的结果精度更高。
(注:
matlab使用的是2015b版本,不知道是保留小数位数太少,还是程序原因,顺序消元输出结果总是等于准确解,请老师指正)
第四次上机作业
7.分析用下列迭代法解线性方程组
=
的收敛性,并求出使
0.0001的近似解及相应的迭代次数。
(1)雅可比迭代法
解:
上机编写的函数如下
function[]=Jacobi(X,b)
%雅可比迭代法
D=diag(diag(X));%得到对角线元素组成的矩阵
B=inv(D)*(D-X);
f=inv(D)*b;
b(:
:
)=0;
x1=B*b+f;
num=1;
while(norm(x1-b)>0.0001)%判断当前的解是否达到精度要求
b=x1;
x1=B*b+f;
num=num+1;
end;
fprintf('求得的解为:
\n');
disp(x1);
fprintf('迭代次数:
%d次\n',num);
end
上机运行,结果如下
>>A=[4,-1,0,-1,0,0;-1,4,-1,0,-1,0;0,-1,4,-1,0,-1;-1,0,-1,4,-1,0;0,-1,0,-1,4,-1;0,0,-1,0,-1,4];
>>b=[0;5;-2;5;-2;6];
>>Jacobi(A,b)
求得的解为:
0.999981765074381
1.99995018125674
0.999975090628368
1.99995018125674
0.999975090628368
1.99996353014876
迭代次数:
28次
满足要求的解如输出结果所示,总共迭代了28次
(2)高斯-赛德尔迭代法
上机程序如下所示
function[]=Gauss_Seidel(A,b)
%高斯赛德尔迭代法
D=diag(diag(A));
L=D-tril(A);
U=D-triu(A);
B=inv(D-L)*U;
f=inv(D-L)*b;
b(:
:
)=0;
x0=B*b+f;
num=1;
while(norm(x0-b)>0.0001)
num=num+1;
b=x0;
x0=B*b+f;
end;
fprintf('结果为\n');
disp(x0);
fprintf('迭代次数为:
%d次\n',num);
end
>>A=[4,-1,0,-1,0,0;-1,4,-1,0,-1,0;0,-1,4,-1,0,-1;-1,0,-1,4,-1,0;0,-1,0,-1,4,-1;0,0,-1,0,-1,4];
>>b=[0;5;-2;5;-2;6];
>>Gauss_Seidel(A,b)
结果为
0.999960143810658
1.99995676152139
0.999963508299833
1.99996662162874
0.999972527179715
1.99998400886989
迭代次数为:
15次
满足精度要求的解如上述程序打印输出所示,迭代了15次
(3)SOR迭代法(w依次取1.334,1.95,0.95)
上机代码如下
function[]=SOR(A,b,w)
%SOR迭代法¨
D=diag(diag(A));
L=D-tril(A);
U=D-triu(A);
B=inv(D-w*L)*((1-w)*D+w*U);
f=w*inv(D-w*L)*b;
b(:
:
)=0;
x0=B*b+f;
num=1;
while(norm(x0-b)>0.0001)
num=num+1;
b=x0;
x0=B*b+f;
end;
fprintf('结果为\n');
disp(x0);
fprintf('迭代次数为%d\n',num);
end
上机运行
>>A=[4,-1,0,-1,0,0;-1,4,-1,0,-1,0;0,-1,4,-1,0,-1;-1,0,-1,4,-1,0;0,-1,0,-1,4,-1;0,0,-1,0,-1,4];
>>b=[0;5;-2;5;-2;6];
>>SOR(A,b,1.334)
结果为
1.00001878481009
1.99998698322858
1.00001815013068
2.00000041318053
0.999991858543476
2.0000068413569
迭代次数为13
>>SOR(A,b,1.95)
结果为
0.999984952088107
2.00000960832604
0.999959115182729
2.0000168426006
1.00000443526697
1.99997885113446
迭代次数为241
>>SOR(A,b,0.95)
结果为
0.999961518309351
1.99995706825231
0.999963054838453
1.99996580572033
0.999971141727589
1.9999827300678
迭代次数为17
由以上输出得到w取值不同的情况下,得到的满足精度要求的结果,迭代次数分别如输出所示
第五次上机作业
8.从函数表
x
0.0
0.1
0.2
0.3
0.401
0.5
f(x)
0.39894
0.39695
0.39142
0.38138
0.36812
0.35206
出发,用下列方法计算f(0.15),f(0.31)及f(0.47)的近似值
(1)分段线性插值
(2)分段二次插值
(3)全区间上拉格朗日插值
解:
(1)线性插值
编写函数如下
function[R]=div_line(x0,y0,x)
%线性插值
p=length(x0);
n=length(y0);
m=length(x);
if(p~=n)%x的个数与y的个数不等
error('数据输入有误,请重新输入');
return;
else
fprintf('线性插值\n');
fort=1:
m
z=x(t);
if(z(1)||z>x0(p))
fprintf('x[%d]不在所给自变量范围内,无法进行插值',t);
continue;
else
fori=1:
p-1
if(zbreak;
end;
end;
R(t)=y0(i)*(x(t)-x0(i+1))/(x0(i)-x0(i+1))+y0(i+1)*(x(t)-x0(i))/(x0(i+1)-x0(i));
end;
end;
end;
end
上机运行如下
>>x0=[0.00.10.1950.30.4010.5];
>>y0=[0.398940.396950.391420.381380.368120.35206];
>>x=[0.150.310.47];
>>div_line(x0,y0,x)
线性插值
ans=
0.394040.380070.35693
即结果为f(0.15)
0.39404,f(0.31)
0.38007,f(0.47)
0.35693
(2)分段二次插值
编写的函数如下
function[R]=div2line(x0,y0,x)
%分段二次插值
p=length(x0);
m=length(y0);
n=length(x);
if(p~=m)
error('输入错误,请重新输入数据');
end;
fort=1:
n
if(x(t)(1)||x(t)>x0(p))
fprintf('x[%d]不在所给区间上',t);
continue;
end;
tag=2;
m=abs(x(t)-x0
(1))+abs(x(t)-x0
(2))+abs(x(t)-x0(3));
fori=3:
p-1
sum=abs(x(t)-x0(i-1))+abs(x(t)-x0(i))+abs(x(t)-x0(i+1));
if(summ=sum;
tag=i;
end
end;
fprintf('tag=%d\n',tag);
R(t)=y0(tag-1)*(x(t)-x0(tag))*(x(t)-x0(tag+1))/((x0(tag-1)-x0(tag))*(x0(tag-1)-x0(tag+1)))+y0(tag)*(x(t)-x0(tag-1))*(x(t)-x0(tag+1))/((x0(tag)-x0(tag-1))*(x0(tag)-x0(tag+1)))+y0(tag+1)*(x(t)-x0(tag-1))*(x(t)-x0(tag))/((x0(tag+1)-x0(tag-1))*(x0(tag+1)-x0(tag)));
End
上机运行,执行结果为
>>x0=[0.00.10.1950.30.4010.5];
>>y0=[0.398940.396950.391420.381380.368120.35206];
>>x=[0.150.310.47];
>>div2line(x0,y0,x)
ans=
0.394480.380220.35725
即分段二次插值方法下,
f(0.15)
0.39448,f(0.31)
0.38022,f(0.47)
0.35725
(3)
上机编写的程序如下
function[R]=lagrange(x0,y0,x)
%全区间上拉格朗日插值
p=length(y0);n=length(x0);m=length(x);
%计算函数表和x的长度
ifp~=nerror('数据输入有误,请重新输入');
%若函数表的x与y长度不一致则输入有误
elsefprintf('拉格朗日插值\n');
fort=1:
m
%利用循环计算每个x的插值
s=0.0;
z=x(t);
fork=1:
n
p=1;
fori=1:
n
ifi~=k
p=p*(z-x0(i))/(x0(k)-x0(i));
end
end
s=s+y0(k)*p;
end
%根据拉格朗日插值公式求解y
R(t)=s;
%输出插值结果
end
end
上机运行结果为
>>x0=[0.00.10.1950.30.4010.5];
>>y0=[0.398940.396950.391420.381380.368120.35206];
>>x=[0.150.310.47];
>>lagrange(x0,y0,x)
拉格朗日插值
ans=
0.394470.380220.35722
即分段二次插值方法下,
f(0.15)
0.39447,f(0.31)
0.38022,f(0.47)
0.35722
9.解:
上机程序如下,为方便起见,将所有操作分在四个函数中进行
入口函数
function[]=spline(X,Y,xx,y1_0,y1_18)
%输出自变量所对应的函数值
M=getM(X,Y,y1_0,y1_18);%先得到M
s=xx;
k=length(xx);
fora=1:
k
s(xx(a))=getExp(X,Y,M,xx(a));%计算自变量所在小区间对应曲线的表达式,并根据表达式计算对应的函数值
fprintf('s(%d)=%f\n',xx(a),s(xx(a)));%输出打印函数值
end
end
获取M
function[M]=getM(X,Y,y1_0,y1_1)
%得到M
n=length(X);
a=0*X;b=a;c=a;h=a;f=a;
b=b+2;
h(2:
n)=X(2:
n)-X(1:
n-1);%h
(1)不用
a(2:
n-1)=h(2:
n-1)./(h(2:
n-1)+h(3:
n));
c(2:
n-1)=1-a(2:
n-1);
a(n)=1;c
(1)=1;
Yf(2:
n)=Y(2:
n)-Y(1:
n-1);
f(2:
n-1)=6*(Yf(3:
n)./h(3:
n)-Yf(2:
n-1)./h(2:
n-1))./(h(2:
n-1)+h(3:
n));
f
(1)=6*(Yf
(2)/h
(2)-y1_0)/h
(2);
f(n)=6*(y1_1-Yf(n)/h(n))/h(n);
M=CalM(n,a,b,c,f);%计算M
end
计算M
function[f]=CalM(n,a,b,c,f)
%追赶法求解M
eps=0.1e-15;%防止参数过小,是的计算误差过大
ifabs(b
(1))disp('除数为0,停止计算');
return
else
c
(1)=c
(1)/b
(1);
end
%追赶法:
根据递推算式计算β
fori=2:
n-1
b(i)=b(i)-a(i)*c(i-1);
ifabs(b(i))disp('除数为0,停止计算');
return
else
c(i)=c(i)/b(i);
end
end
b(n)=b(n)-a(n)*c(n-1);
%追赶法:
根据递推算式计算
f
(1)=f
(1)/b
(1);
fori=2:
n
f(i)=(f(i)-a(i)*f(i-1))/b(i);
end
%以下求解Ux=y,x的值存入f
fori=n-1:
-1:
1
f(i)=f(i)-c(i)*f(i+1);
end
return
end
得到自变量所在区间的表达式,并求自变量对应的函数值
function[y]=getExp(X,Y,M,x)
%%根据X、Y、M计算表达式,并根据表达式计算对应的函数值
n=length(X);
h(2:
n)=X(2:
n)-X(1:
n-1);
%%判断x落在哪个小区间
n1=1;n2=n;
whilen2~=n1+1
n5=fix((n1+n2)/2);
ifx>X(n5)
n1=n5;
else
n2=n5;
end
end
%%
%%计算y
y=M(n1)*(X(n2)-x)^3/(6*h(n2))+M(n2)*(x-X(n1))^3/(6*h(n2));
y=y+(Y(n1)-M(n1)*h(n2)*h(n2)/6)*(X(n2)-x)/h(n2);
y=y+(Y(n2)-M(n2)*h(n2)*h(n2)/6)*(x-X(n1))/h(n2);
%%
%结束
end
上机试运行,过程如下
>>X=[0.523.18.017.9528.6539.6250.6578104.6156.6208.6260.7312.5364.4416.3468494507520];
>>Y=[5.287949.413.8420.224.928.4431.13536.536.634.631.026.3420.914.87.83.71.50.2];
>>xx=[24612163060110180280400515];
>>y1_0=1.86548;
>>y1_18=-0.046115;spline(X,Y,xx,y1_0,y1_18)
s
(2)=7.825123
s(4)=10.481311
s(6)=12.363477
s(12)=16.575574
s(16)=19.091594
s(30)=25.386402
s(60)=32.804283
s(110)=36.647878
s(180)=35.917148
s(280)=29.368462
s(400)=16.799784
s(515)=0.542713
根据上述程序运行结果,可得所有自变量对应的函数值,如上输出结果所示
第六次上机作业
10.已知一组实验数据
i
1
2
3
4
5
6