计算方法上机作业集合1Word文件下载.docx

计算方法上机作业集合1Word文件下载.docx

《计算方法上机作业集合1Word文件下载.docx》由会员分享，可在线阅读，更多相关《计算方法上机作业集合1Word文件下载.docx（31页珍藏版）》请在冰豆网上搜索。

上机作业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>

disp（'

此方程组无解'

）;

return

end

ifRA==RB

ifRA==n

formatlong;

此方程组有唯一解'

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

此方程组有无数组解'

上机运行结果为

A=[1.1348,3.8326,1.1651,3.4017;

X=gaus（A,b）

此方程组有唯一解

X=

1.000000000000000

（2）列主元消元法

函数部分代码如下

function[b]=lieZhu（A,b）%用b返回方程组的解

formatlong;

该方程组无解'

ifdif==0

该方程组有唯一解'

c=zeros（1,n）;

fori=1:

max=abs（A（i,i））;

m=i;

forj=i+1:

ifmax<

abs（A（j,i））

max=abs（A（j,i））;

m=j;

end%求出每一次消元时绝对值最大的一行的行号

ifm~=i

fork=i:

c（k）=A（i,k）;

A（i,k）=A（m,k）;

A（m,k）=c（k）;

d1=b（i）;

b（i）=b（m）;

b（m）=d1;

%函数值跟随方程一起换位置

end

fork=i+1:

A（k,j）=A（k,j）-A（i,j）*A（k,i）/A（i,i）;

b（k）=b（k）-b（i）*A（k,i）/A（i,i）;

A（k,i）=0;

end%完成消元操作，形成上三角矩阵

fori=n-1:

sum=0;

sum=sum+A（i,j）*b（j）;

b（i）=（b（i）-sum）/A（i,i）;

%回代求解其他未知数

else

end上机运行，结果为

X=lieZhu（A,b）

该方程组有唯一解

1.000000000000002

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

迭代次数：

%d次\n'

num）;

上机运行，结果如下

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;

6];

Jacobi（A,b）

0.999981765074381

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;

x0=B*b+f;

while（norm（x0-b）>

0.0001）

b=x0;

x0=B*b+f;

结果为\n'

disp（x0）;

迭代次数为：

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迭代法¨

B=inv（D-w*L）*（（1-w）*D+w*U）;

f=w*inv（D-w*L）*b;

迭代次数为%d\n'

上机运行

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;

fprintf（'

线性插值\n'

fort=1:

m

z=x（t）;

if（z<

x0

（1）||z>

x0（p））

x[%d]不在所给自变量范围内，无法进行插值'

t）;

continue;

p-1

x0（i+1））

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

上机运行如下

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）

%分段二次插值

m=length（y0）;

n=length（x）;

if（p~=m）

输入错误，请重新输入数据'

fort=1:

if（x（t）<

x0

（1）||x（t）>

x[%d]不在所给区间上'

tag=2;

m=abs（x（t）-x0

（1））+abs（x（t）-x0

（2））+abs（x（t）-x0（3））;

fori=3:

sum=abs（x（t）-x0（i-1））+abs（x（t）-x0（i））+abs（x（t）-x0（i+1））;

if（sum<

m）

m=sum;

tag=i;

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

上机运行，执行结果为

div2line（x0,y0,x）

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

%计算函数表和x的长度

ifp~=nerror（'

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

elsefprintf（'

拉格朗日插值\n'

m

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

s=0.0;

fork=1:

n

p=1;

ifi~=k

p=p*（z-x0（i））/（x0（k）-x0（i））;

s=s+y0（k）*p;

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

R（t）=s;

%输出插值结果

lagrange（x0,y0,x）

拉格朗日插值

0.394470.380220.35722

0.39447,f（0.31）

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

%计算自变量所在小区间对应曲线的表达式，并根据表达式计算对应的函数值

s（%d）=%f\n'

xx（a）,s（xx（a）））;

%输出打印函数值

获取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:

a（n）=1;

c

（1）=1;

Yf（2:

n）=Y（2:

n）-Y（1:

f（2:

n-1）=6*（Yf（3:

n）./h（3:

n）-Yf（2:

n-1）./h（2:

n-1））./（h（2:

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

计算M

function[f]=CalM（n,a,b,c,f）

%追赶法求解M

eps=0.1e-15;

%防止参数过小，是的计算误差过大

ifabs（b

（1））<

eps

除数为0，停止计算'

c

（1）=c

（1）/b

（1）;

%追赶法：

根据递推算式计算β

fori=2:

b（i）=b（i）-a（i）*c（i-1）;

ifabs（b（i））<

c（i）=c（i）/b（i）;

b（n）=b（n）-a（n）*c（n-1）;

%追赶法:

根据递推算式计算

f

（1）=f

（1）/b

（1）;

f（i）=（f（i）-a（i）*f（i-1））/b（i）;

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

fori=n-1:

f（i）=f（i）-c（i）*f（i+1）;

return

得到自变量所在区间的表达式，并求自变量对应的函数值

function[y]=getExp（X,Y,M,x）

%%根据X、Y、M计算表达式，并根据表达式计算对应的函数值

%%判断x落在哪个小区间

n1=1;

n2=n;

whilen2~=n1+1

n5=fix（（n1+n2）/2）;

ifx>

X（n5）

n1=n5;

n2=n5;

%%

%%计算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）;

%结束

上机试运行，过程如下

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

2

3

4

5

6