计算方法上机作业集合1.docx

计算方法上机作业集合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

ifmax

max=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（z

break;

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