计算方法上机报告.docx

计算方法上机报告.docx

《计算方法上机报告.docx》由会员分享，可在线阅读，更多相关《计算方法上机报告.docx（28页珍藏版）》请在冰点文库上搜索。

计算方法上机报告

数值计算方法上机实习

报告

1.设

，

（1）由递推公式

，从

的几个近似值出发，计算

；

（2）粗糙估计

，用

，计算

；

（3）分析结果的可靠性及产生此现象的原因（重点分析原因）。

解：

（1）易得：

log（6）-log（5）=0.1823,

程序为：

I=0.1823;

forn=1:

20

I=（-5）*I+1/n;

end

I

输出结果为：

=-2.0558e+009

（2）粗糙估计

，用

，计算

；

=（1/6/21+1/5/21）/2=0.0087,

程序为：

I=0.0087;

forn=20:

-1:

1

I=（-1/5）*I+1/（5*n）;

end

I

输出结果：

I=0.1823

（3）分析结果的可靠性及产生此现象的原因（重点分析原因）。

首先分析两种递推式的误差；设第一递推式中开始时的误差为

，递推过程的舍入误差不计。

并记

，则有

。

因为

，所此递推式不可靠。

而在第二种递推式中

，误差在缩小，所以此递推式是可靠的。

出现以上运行结果的主要原因是在构造递推式过程中，考虑误差是否得到控制，即算法是否数值稳定。

2.求方程

的近似根，要求

，并比较计算量。

（1）在[0,1]上用二分法；

（2）取初值

，并用迭代

；

（3）加速迭代的结果；

（4）取初值

，并用牛顿迭代法；

（5）分析绝对误差。

（1）在[0，1]上用二分法；

程序：

a=0;b=1.0;

whileabs（b-a）>5*1e-4

c=（b+a）/2;

ifexp（c）+10*c-2>0

b=c;

elsea=c;

end

end

c

结果：

c=0.0903

（2）取初值

，并用迭代

；

程序：

x=0;

a=1;

whileabs（x-a）>5*1e-4

a=x;

x=（2-exp（x））/10;

end

x

结果：

x=0.0905

（3）加速迭代的结果；

程序：

x=0;

a=0;b=1;

whileabs（b-a）>5*1e-4

a=x;

y=（2-exp（x））/10;

z=（2-exp（y））/10;

x=x-（y-x）^2/（z-2*y+x）;

b=x;

end

x

结果：

x=0.0905

（4）取初值

，并用牛顿迭代法；

程序：

x=0;

a=0;b=1;

whileabs（b-a）>5*1e-4

a=x;

x=x-（exp（x）+10*x-2）/（exp（x）+10）;

b=x;

end

x

结果：

x=0.0905

（5）分析绝对误差。

x=solve（'exp（x）+10*x-2=0'）

可得：

x=0.090525

可知二分法绝对误差最大，不动点法迭代、加速迭代和牛顿法效果都好。

3.钢水包使用次数多以后，钢包的容积增大，数据如下：

x

2

3

4

5

6

7

8

9

y

6.42

8.2

9.58

9.5

9.7

10

9.93

9.99

10

11

12

13

14

15

16

10.49

10.59

10.60

10.8

10.6

10.9

10.76

试从中找出使用次数和容积之间的关系，计算均方差。

（注：

增速减少，用何种模型）

解：

设y=f（x）具有指数形式

（a>0,b<0）。

对此式两边取对数，得

。

记A=lna，B=b，并引入新变量z=lny，t=1/x。

引入新变量后的数据表如下

x

2

3

4

5

6

7

8

9

t=1/x

0.5000

0.3333

0.2500

0.2000

0.1667

0.1429

0.1250

0.1111

z=lny

1.8594

2.1041

2.2597

2.2513

2.2721

2.3026

2.2956

2.3016

10

11

12

13

14

15

16

0.1000

0.0909

0.0833

0.0769

0.0714

0.0667

0.0625

2.3504

2.3599

2.3609

2.3795

2.3609

2.3888

2.3758

程序：

t=[0.50000.33330.25000.20000.16670.14290.12500.11110.10000.09090.08330.07690.07140.06670.0625];

z=[1.85942.10412.25972.25132.27212.30262.29562.30162.35042.35992.36092.37952.36092.38882.3758];

polyfit（t,z,1）

结果：

ans=-1.11072.4578

由此可得A=2.4578，B=-1.1107，

，b=B=-1.1107

方程即为

计算均方差编程：

x=[2:

16];y=[6.428.29.589.59.7109.939.9910.4910.5910.6010.810.610.910.76];

f（x）=11.6791*exp（-1.1107./x）;

c=0;

fori=1:

15

a=y（i）;

b=x（i）;

c=c+（a-f（b））^2;

end

averge=c/15

结果：

averge=0.0594

4.设

，

，

分析下列迭代法的收敛性，并求

的近似解及相应的迭代次数。

（1）JACOBI迭代；

（2）GAUSS-SEIDEL迭代；

（3）SOR迭代（

）。

（1）解：

JACOBI迭代；

程序如下：

function[x,k,index]=Jacobi（A,b,ep,it_max）

%求解线性方程组的Jacobi迭代法，其中

%A---方程组的系数矩阵

%b---方程组的右端项

%ep---精度要求。

省缺为1e-4

%it_max---最大迭代次数，省缺为100

%x---方程组的解

%k---迭代次数

%index---index=1表示迭代收敛到指定要求；

%index=0表示迭代失败

n=length（A）;k=0;

x=zeros（n,1）;y=zeros（n,1）;index=1;

while1

fori=1:

n

y（i）=b（i）;

forj=1:

n

ifj~=i

y（i）=y（i）-A（i,j）*x（j）;

end

end

ifk==it_max

index=0;return;

end

y（i）=y（i）/A（i,i）;

end

ifnorm（y-x）

break;

end

x=y;k=k+1;

end

文件以Jacobi.m保存

程序：

A=[4-10-100;-14-10-10;0-14-10-1;-10-14-10;0-10-14-1;00-10-14];

b=[05-25-26]';

[x,k,index]=Jacobi（A,b,1e-4,100）

运行结果为：

x=

0.999946604431898

1.99996353014876

0.999927060297523

1.99996353014876

0.999927060297523

1.99997330221595

k=

27

index=

1

（2）GAUSS-SEIDEL迭代；

function[x,k,flag]=gau_seid（A,b,delta,maxl）

%求解线性方程组的迭代法，其中

%A为方程组的系数矩阵

%b为方程组的右端项

%delta为代数精度，默认值为1e-4

%maxl为最大迭代次数，默认值为100

%x为方程组的解；

%k为迭代次数

%flag为指标变量flag='OK!

'表示迭代收敛指标要求

%flag='fail'表示迭代失败

n=length（A）;

k=0;

x=zeros（n,1）;

y=zeros（n,1）;

flag='OK!

';

whilek<100

y=x;

fori=1:

n

z=b（i）;

forj=1:

n

ifj~=i

z=z-A（i,j）*x（j）;

end

end

%ifk==maxl

%flag='Fail';

%return

%end

z=z/A（i,i）;

x（i）=z;

end

ifnorm（y-x）

break

end

k=k+1;

end

以文件名gau_seid.m保存。

程序：

A=[4-10-100;-14-10-10;0-14-10-1;-10-14-10;0-10-14-1;00-10-14];

b=[05-25-26]';

[x,k,flag]=gau_seid（A,b,1e-4,100）

运行结果为：

x=

0.999960143810659

1.99995676152139

0.999963508299833

1.99996662162874

0.999972527179715

1.99998400886989

k=

14

flag=

OK!

（3）SOR迭代（

）。

程序：

function[x,k,flag]=sor（A,b,delta,w,maxl）

%求解线性方程组的迭代法，其中

%A为方程组的系数矩阵

%b为方程组的右端项

%delta为代数精度，默认值为1e-4

%w为松弛因子

%maxl为最大迭代次数，默认值为100

%x为方程组的解；

%k为迭代次数

%flag为指标变量flag='OK!

'表示迭代收敛指标要求

%flag='fail'表示迭代失败

n=length（A）;

k=0;

x=zeros（n,1）;

y=zeros（n,1）;

flag='OK!

';

whilek<100

y=x;

fori=1:

n

z=b（i）;

forj=1:

n

ifj~=i

z=z-A（i,j）*x（j）;

end

end

%ifk==maxl

%flag='Fail';

%return

%end

z=z/A（i,i）;

x（i）=（1-w）*x（i）+w*z;

end

ifnorm（y-x）

break

end

k=k+1;

end

以文件名sor.m保存。

程序：

A=[4-10-100;-14-10-10;0-14-10-1;-10-14-10;0-10-14-1;00-10-14];

b=[05-25-26]';

[x,k,flag]=sor（A,b,1e-4,1.334,100）

[x,k,flag]=sor（A,b,1e-4,1.95,100）

[x,k,flag]=sor（A,b,1e-4,0.95,100）

运行结果分别为：

x=

1.00001878481009

1.99998698322858

1.00001815013068

2.00000041318053

0.999991858543476

2.0000068413569

k=

12

flag=

OK!

x=

1.00708190171991

2.00285345528501

1.01400422955242

1.98988871374656

1.00950321530473

2.00707099515554

k=

100

flag=

OK!

x=

0.999961518309351

1.99995706825231

0.999963054838452

1.99996580572033

0.999971141727588

1.9999827300678

k=

16

flag=

OK!

5.用逆幂迭代法求

最接近于11的特征值和特征向量，准确到

。

解：

算法为结合原点平移的反幂法，编写MATLAB程序如下：

functiont=maxnorm（A）

n=length（A）;

t=0;

fori=1:

n

ifabs（A（i））-abs（t）>=0

t=A（i）;

end

end

function函数保存为maxnorm.m文件

程序：

A=[631;321;111];

v=[111]';

lam

（1）=maxnorm（v）;

u=v/lam

（1）;

fori=1:

20

fprintf（'n=%d,x=[%f%f%f],lamda=%f\n',i-1,u

（1）,u

（2）,u（3）,1/lam（i）+11）;

v=（A-11*eye（3））\u;

lam（i+1）=maxnorm（v）;

u=v/lam（i+1）;

ifabs（1/lam（i+1）-1/lam（i））<1e-3;

break

end

end

运行结果如下：

n=0.000000,x=[1.0000001.0000001.000000],lamda=12.000000

n=1.000000,x=[1.0000000.6666670.424242],lamda=8.424242

n=2.000000,x=[1.0000000.5843680.284130],lamda=8.037233

n=3.000000,x=[1.0000000.5601180.243417],lamda=7.923772

n=4.000000,x=[1.0000000.5526220.230993],lamda=7.888859

n=5.000000,x=[1.0000000.5502720.227144],lamda=7.877961

n=6.000000,x=[1.0000000.5495330.225946],lamda=7.874545

n=7.000000,x=[1.0000000.5493000.225573],lamda=7.873473

故所求特征值为7.873473，特征向量为[1.0000000.5493000.225573]。

6.用经典R-K方法求解初值问题

（1）

，

，

；

（2）

，

，

。

和精确解

比较，分析结论。

解：

（1）编写MATLAB程序如下：

四阶R-K方法：

function[x,y,z]=nark4（dyfun1,dyfun2,xspan,y0,z0,h）

%y'=f（x,y,z）,z'=g（x,y,z）,y（x0）=y0,z（x0）=z0

%dyfun为函数f（x,y）,xspan为求解区间[x0,xN],y0为初值y（x0）,z0为初值z（x0）,h为步长

%x返回节点，y返回数值解

x=xspan

（1）:

h:

xspan

（2）;

y

（1）=y0;

z

（1）=z0;

forn=1:

length（x）-1

k1=feval（dyfun1,x（n）,y（n）,z（n））;

l1=feval（dyfun2,x（n）,y（n）,z（n））;

k2=feval（dyfun1,x（n）+h/2,y（n）+h/2*k1,z（n）+h/2*l1）;

l2=feval（dyfun2,x（n）+h/2,y（n）+h/2*k1,z（n）+h/2*l1）;

k3=feval（dyfun1,x（n）+h/2,y（n）+h/2*k2,z（n）+h/2*l2）;

l3=feval（dyfun2,x（n）+h/2,y（n）+h/2*k2,z（n）+h/2*l2）;

k4=feval（dyfun1,x（n+1）,y（n）+h*k3,z（n）+h*l3）;

l4=feval（dyfun2,x（n+1）,y（n）+h*k3,z（n）+h*l3）;

y（n+1）=y（n）+h*（k1+2*k2+2*k3+k4）/6;

z（n+1）=z（n）+h*（l1+2*l2+2*l3+l4）/6;

end

程序保存为nark4.m文件。

程序：

dyfun1=inline（'-2*y+z+2*sin（x）'）;

dyfun2=inline（'y-2*z+2*cos（x）-2*sin（x）'）;

[x,y,z]=nark4（dyfun1,dyfun2,[0,10],2,3,0.01）;

u=2*exp（-x）+sin（x）;%u为y的真值

v=2*exp（-x）+cos（x）;%v为z的真值

fori=1:

length（x）

err1=（y（i）-u（i））;

err2=（z（i）-v（i））;

plot（x（i）,err1,'r-'）

holdon

plot（x（i）,err2,'g-'）

end

运行结果：

（2）二阶R-K方法：

function[x,y,z]=nark2（dyfun1,dyfun2,xspan,y0,z0,h）

%y'=f（x,y,z）,z'=g（x,y,z）,y（x0）=y0,z（x0）=z0

%dyfun为函数f（x,y）,xspan为求解区间[x0,xN],y0为初值y（x0）,z0为初值z（x0）,h为步长

%x返回节点，y返回数值解

x=xspan

（1）:

h:

xspan

（2）;

y

（1）=y0;

z

（1）=z0;

forn=1:

length（x）-1

k1=feval（dyfun1,x（n）,y（n）,z（n））;

I1=feval（dyfun2,x（n）,y（n）,z（n））;

k2=feval（dyfun1,x（n）+h,y（n）+h*k1,z（n）+h*I1）;

I2=feval（dyfun2,x（n）+h,y（n）+h*k1,z（n）+h*I1）;

y（n+1）=y（n）+h*（k1/2+k2/2）;

z（n+1）=z（n）+h*（I1/2+I2/2）;

end

程序保存为nark2.m文件。

程序：

dyfun1=inline（'-2*y+z+2*sin（x）'）;

dyfun2=inline（'998*y-999*z+999*cos（x）-999*sin（x）'）;

[x,y,z]=nark2（dyfun1,dyfun2,[0,10],2,3,0.0005）;

u=2*exp（-x）+sin（x）;%u为y的真值

v=2*exp（-x）+cos（x）;%v为z的真值

fori=1:

length（x）

err1=（y（i）-u（i））;

err2=（z（i）-v（i））;

plot（x（i）,err1,'r-'）

holdon

plot（x（i）,err2,'g-'）

end

运行结果：

7.用有限差分法求解边值问题（h=0.1）：

。

解：

通过差商逼近导数，将微分方程离散为差分方程

编写MATLAB程序如下：

a0=-1;b0=1;c0=1;d0=1;

h=0.1;

fori=1:

19;

x=a0+i*h;

q=-（1+x*x）;

a（i）=1;b（i）=-2+h*h*q;c（i）=1;d（i）=0;

fprintf（'%f%f%f%f\n',a（i）,b（i）,c（i）,d（i））

end

运行结果如下：

1.000000-2.0181001.0000000.000000

1.000000-2.0164001.0000000.000000

1.000000-2.0149001.0000000.000000

1.000000-2.0136001.0000000.000000

1.000000-2.0125001.0000000.000000

1.000000-2.0116001.0000000.000000

1.000000-2.0109001.0000000.000000

1.000000-2.0104001.0000000.000000

1.000000-2.0101001.0000000.000000

1.000000-2.0100001.0000000.000000

1.000000-2.0101001.0000000.000000

1.000000-2.0104001.0000000.000000

1.000000-2.0109001.0000000.000000

1.000000-2.0116001.0000000.000000

1.000000-2.0125001.0000000.000000

1.000000-2.0136001.0000000.000000

1.000000-2.0149001.0000000.000000

1.000000-2.0164001.0000000.000000

1.000000-2.0181001.0000000.000000

采用追赶法解该线性方程组，编写MATLAB程序如下：

a=[1111111111111111111];

b=[-2.0181-2.0164-2.0149-2.0136-2.0125-2.0116-2.0109-2.0104-2.0101-2.0100-2.0101-2.0104-2.0109-2.0116-2.0125-2.0136-2.0149-2.0164-2.0181];

c=[1111111111111111111];

d=[-100000000000000000-1];

x=d;

n=length（x）;

y

（1）=1;

y（21）=1;

forj=1:

n-1

mu=a（j）/b（j）;

b（j+1）=b（j+1）-mu*c（j）;

x（j+1）=x（j+1）-mu*x（j）;

end

x（n）=x（n）/b（n）;

forj=n-1:

-1:

1

x（j）=（x（j）-c（j）*x（j+1））/b（j）;

end

form=1:

1:

19

y（m+1）=x（m）;

end

h=-1:

0.1:

1;

plot（h,y）;

grid