计算方法上机作业集合

第一次&第二次上机作业

上机作业：

1.在Matlab上执行：>> 5.1-5-0.1和>> 1.5-1-0.5

给出执行结果，并简要分析一下产生现象的原因。

解：执行结果如下：

在Matlab中，小数值很难用二进制进行描述。由于计算精度的影响，相近两数相减会出现误差。

2.（课本181页第一题）

解：（1）n=0时，积分得I0=ln6-ln5，编写如下图代码

从以上代码显示的结果可以看出，I 20的近似值为0.7465

（2）I I =∫I I 5+I 10dx,可得∫I I 610dx ≤∫I I 5+I 10dx ≤∫I I

510dx,得 16(I +1)≤I I ≤15(I +1),则有1126≤I 20≤1105,

取I 20=1

105

,以此逆序估算I 0。代码段及结果如下图所示

（3）从I20估计的过程更为可靠。首先根据积分得表达式是可知，被积函数随着n的增大，其所围面积应当是逐步减小的，即积分值应是随着n的递增二单调减小的，（1）中输出的值不满足这一条件，（2）满足。设I I表示I I的近似值，I I-I I=(?5)I(I0?I0)(根据递推公式可以导出此式),可以看出，随着n的增大，误差也在增大，所以顺序估计时，算法不稳定性逐渐增大，逆序估计情况则刚好相反，误差不断减小，算法逐渐趋于稳定。

2.（课本181页第二题）

（1）上机代码如图所示

求得近似根为0.09058 （2）上机代码如图所示

得近似根为0.09064；

（3）牛顿法上机代码如下

计算所得近似解为0.09091

第三次上机作业上机作业181页第四题

线性方程组为

[1.13483.8326

0.53011.7875

1.16513.4017

2.53301.5435

3.4129

4.9317

1.23714.9998

8.76431.3142

10.67210.0147

][

]=[

9.5342

6.3941

18.4231

16.9237

]

（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;

if dif>0

disp('此方程组无解');

return

end

if RA==RB

if RA==n

format long;

disp('此方程组有唯一解');

for p=1:n-1

for k=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);

for q=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

上机运行结果为

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.0000

（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;

format long;

if dif>0

disp('该方程组无解');

return

end

if dif==0

if RA==n

disp('该方程组有唯一解');

c=zeros(1,n);

for i=1:n-1

max=abs(A(i,i));

m=i;

for j=i+1:n

if max

max=abs(A(j,i));

m=j;

end

end %求出每一次消元时绝对值最大的一行的行号 if m~=i

for k=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

for k=i+1:n

for j=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);

for i=n-1:-1:1

sum=0;

for j=i+1:n

sum=sum+A(i,j)*b(j);

end

b(i)=(b(i)-sum)/A(i,i) ;%回代求解其他未知数

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.0000

1.0002

0.9999

根据两种方法运算结果，顺序消元法得到结果比列主元消元法得到的结果精度更高。

（注：matlab使用的是2015b版本，不知道是保留小数位数太少，还是程序原因，顺序消元输出结果总是等于准确解，请老师指正）

第四次上机作业

7.分析用下列迭代法解线性方程组

[

4?1

?14

0?1

?10

0?1

?10

4?1

?14

0 ?1

?10

0?1

?10

4?1

?14][

I6]

[

的收敛性，并求出使‖I(I+1)?I(I)‖

≤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.4381

1.674

0.8368

1.674

0.8368

1.876

迭代次数：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.0658

1.139

0.9833

1.874

0.9715

1.989

迭代次数为：15次

满足精度要求的解如上述程序打印输出所示，迭代了15次（3）S OR迭代法（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.009

1.858

1.068

2.053

0.3476

2.69

迭代次数为13

>> SOR(A,b,1.95)

结果为

0.8107

2.604

0.2729

2.06

1.697

1.446

迭代次数为241

>> SOR(A,b,0.95)

结果为

0.9351

1.231

0.8453

1.033

0.7589

1.78

迭代次数为17

由以上输出得到w取值不同的情况下，得到的满足精度要求的结果，迭代次数分别如输出所示

第五次上机作业

8.从函数表

出发，用下列方法计算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');

for t=1:m

z=x(t);

if(zx0(p))

fprintf('x[%d]不在所给自变量围，无法进行插值',t);

continue;

else

for i=1:p-1

if(z

break;

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

上机运行如下

>> x0=[0.0 0.1 0.195 0.3 0.401 0.5];

>> y0=[0.39894 0.39695 0.39142 0.38138 0.36812

0.35206];

>> x=[0.15 0.31 0.47];

>> div_line(x0,y0,x)

线性插值

ans =

0.39404 0.38007 0.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;

for t=1:n

if(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));

for i=3:p-1

sum=abs(x(t)-x0(i-1))+abs(x(t)-x0(i))+abs(x(t)-x0(i+1));

if(sum

m=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(t

ag+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.0 0.1 0.195 0.3 0.401 0.5];

>> y0=[0.39894 0.39695 0.39142 0.38138 0.36812 0.35206];

>> x=[0.15 0.31 0.47];

>>div2line(x0,y0,x)

ans =

0.39448 0.38022 0.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的长度

if p ~= n error('数据输入有误，请重新输入');

%若函数表的x与y长度不一致则输入有误

else fprintf('拉格朗日插值\n');

for t=1:m

%利用循环计算每个x的插值

s=0.0;

z=x(t);

for k=1:n

p=1;

for i=1:n

if i~=k

p=p*(z-x0(i))/(x0(k)-x0(i));

end

s=s+y0(k)*p;

end

%根据拉格朗日插值公式求解y

R(t)=s;

%输出插值结果

end

上机运行结果为

>> x0=[0.0 0.1 0.195 0.3 0.401 0.5];

>> y0=[0.39894 0.39695 0.39142 0.38138 0.36812

0.35206];

>> x=[0.15 0.31 0.47];

>> lagrange(x0,y0,x)

拉格朗日插值

ans =

0.39447 0.38022 0.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);

for a=1:k

s(xx(a))=getExp(X,Y,M,xx(a));%计算自变量所在小区间对应曲线的表达式，并根据表达式计算对应的函数值

fprintf('s(%d)=%f\n',xx(a),s(xx(a))); %输出打印函数值

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; %防止参数过小，是的计算误差过大

if abs(b(1))

disp('除数为0，停止计算');

return

else

c(1)=c(1)/b(1);

end

%追赶法：根据递推算式计算β

for i=2:n-1

b(i)=b(i)-a(i)*c(i-1);

if abs(b(i))

disp('除数为0，停止计算');

return

else

c(i)=c(i)/b(i);

end

b(n)=b(n)-a(n)*c(n-1);

%追赶法:根据递推算式计算

f(1)=f(1)/b(1);

for i=2:n

f(i)=(f(i)-a(i)*f(i-1))/b(i);

end

%以下求解Ux=y, x的值存入f

for i=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;

while n2~=n1+1

n5=fix((n1+n2)/2);

if x>X(n5)

n1=n5;

else

n2=n5;

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);