数值计算方法上机实习题.docx

数值计算方法上机实习题.docx

《数值计算方法上机实习题.docx》由会员分享，可在线阅读，更多相关《数值计算方法上机实习题.docx（21页珍藏版）》请在冰豆网上搜索。

数值计算方法上机实习题

数值计算方法上机实习题

1．设

，

（1）由递推公式

，从

出发，计算

；

（2）

，

用

，计算

；

（1）由I0计算I20

递推子程序：

functionf=fib（n,i）

ifn>=1

f=fib（n-1,i）*（-5）+（1/（n））;

elseifn==0

f=i;

end

计算和显示程序：

I=0.1824;

I1=0.1823;

fib1=fib（20,I）;

fib2=fib（20,I1）;

fprintf（'I_0=0.1824时，I_20=%d\n',fib1）;

fprintf（'I_0=0.1823时，I_20=%d\n',fib2）;

计算结果显示：

I_0=0.1824时，I_20=7.480927e+09

I_0=0.1823时，I_20=-2.055816e+09

（2）由I20计算I0

程序：

n=21;

i1=0;

i2=10000;

f1=i1;

f2=i2;

whilen~=0;

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

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

n=n-1;

end

fprintf（'I_20=0时，I_0=%4.8f\n',f1）;

fprintf（'I_20=10000时，I_0=%4.8f\n',f2）;

计算结果显示：

I_20=0时，I_0=0.18232156

I_20=10000时，I_0=0.18232156

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

答：

第一个算法可得出

易知第一个算法每一步计算都把误差放大了5倍，n次计算后更是放大了5n倍，可靠性低。

第二个算法可得出

可以看出第二个算法每一步计算就把误差缩小5倍，n次后缩小了5n倍，可靠性高。

2．求方程

的近似根，要求

，并比较计算量。

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

（1）[0，1]上的二分法

二分法子程序：

function[root,n]=EFF3（f,x1,x2）

%第二题

（1）二分法

f1=subs（f,symvar（f）,x1）;%函数在x=x1的值

f2=subs（f,symvar（f）,x2）;%x=x2

n=0;%步数

er=5*10^-4;%误差

if（f1==0）

root=x1;

return;

elseif（f2==0）

root=x2;

return;

elseif（f1*f2>0）

disp（'两端点函数值乘积大于0！

'）;

return;

else

while（abs（x1-x2）>er）%循环

x3=（x1+x2）/2;

f3=subs（f,symvar（f）,x3）;

n=n+1;

if（f3==0）

root=x3;

break;

elseif（f1*f3>0）

x1=x3;

else

x2=x3;

end

end

root=（x1+x2）/2;%while循环少一步需加上

end

计算根与步数程序：

fplot（@（x）exp（x）+10*x-2,[0,1]）;

gridon;

symsx;

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

[root,n]=EFF3（f,0,1）;

fprintf（'root=%6.8f,n=%d\n',root,n）;

计算结果显示：

root=0.09057617,n=11

（2）取初值

，并用迭代

；

（2）初值x0=0迭代

迭代法子程序：

function[root,n]=DDF（g,x0,err,max）（接下页）

%root根，n+1步数，g函数，x0初值，err误差，max最大迭代次数

X

（1）=x0;

forn=2:

max

X（n）=subs（g,symvar（g）,X（n-1））;

c=abs（X（n）-X（n-1））;

root=X（n）;

if（c

break;

end

ifn==max

disp（'超出预设迭代次数'）;

end

end

n=n-1;%步数减1

计算根与步数程序：

symsx;

f=（2-exp（x））/10;（接下页）

x0=0;

err=5*10^（-4）;

max=100;

[root,n]=DDF（f,x0,err,max）;

fprintf（'root=%6.8f,n=%d\n',root,n）;

计算结果显示：

root=0.09051262,n=4

（3）加速迭代的结果;

（3）加速迭代

加速迭代计算程序：

x0=0;

err=5*10^（-4）;

max=100;

symsx;

g=x-（f（x）-x）^2/（f（f（x））-2*f（x）-x）;

[root,n]=DDF（g,x0,err,max）;

fprintf（'root=%6.8f,n=%d',root,n）;

程序函数设置：

functiong=f（x）

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

计算结果显示：

root=0.09048375,n=2

（4）取初值

，并用牛顿迭代法；

（4）牛顿迭代法

牛顿迭代法子程序：

function[root,n]=NDDDFx（g,x0,err,max）

%root根，n+1步数，g函数，x0初值，err误差，max最大迭代次数

X

（1）=x0;

forn=2:

max

X（n）=subs（g,symvar（g）,X（n-1））;

c=abs（X（n）-X（n-1））;

root=X（n）;

if（c

break;

end

ifn==max

disp（'超出预设迭代次数'）;（接下页）

end

end

n=n-1;%步数减1

牛顿迭代法计算程序：

x0=0;

err=5*10^（-4）;

max=100;

symsx;

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

g=x-（g/diff（g））;

[root,n]=NDDDFx（g,x0,err,max）;

fprintf（'root=%6.8f,n=%d\n',root,n）；

计算结果显示：

root=0.09052511,n=2

（5）分析绝对误差。

答：

可以看到，在同一精度下，二分法运算了11次，题设迭代算式下运算了4次，加速迭代下运算了2次，牛顿迭代下运算了2次。

因不动点迭代法和二分法都是线性收敛的，但二分法压缩系数比题设迭代方法大，收敛速度较慢。

加速迭代速度是超线性收敛，牛顿法是二阶，收敛速度快。

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

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

（用

拟合）

拟合曲线程序：

x=2:

16;

y=[6.428.29.589.59.7109.939.9910.4910.5910.6010.810.610.910.76];

[f,fval]=fminsearch（@delta1,[0,0,0],optimset,x,y）;

%fminsearch为求解多元函数的最小值函数

%f为多元函数初值x0附近的极小值点

%fval为极小值

k=2:

100;

y1=（f

（1）.*k+f

（2））./（f（3）+k）;

figure1=figure（'Color',[111]）;

gxt=plot（x,y,'xr',k,y1,'k-'）;

legend（'原数据','拟合曲线','Location','northwest'）;%y为数据点连接曲线，y1为拟合曲线

title（'函数y=（ax+b）/（x+c）的拟合','FontSize',14,'FontWeight','Bold'）;

xlabel（'次数','FontSize',14,'FontWeight','Bold'）;

ylabel（'容积','FontSize',14,'FontWeight','Bold'）;

set（gxt,'LineWidth',1.5）;

gridon;

%计算均方差

fori=1:

15

y2（i）=（f

（1）.*x（i）+f

（2））./（f（3）+x（i））;（接下页）

end

j=0;

fori=1:

15

j=j+（y（i）-y2（i））^2;

end

jfc=sqrt（j/15）;

fprintf（'拟合出的方程为:

（x+%4.4f）y=%4.4fx+%4.4f\n均方差为:

%4.8f\n',f（3）,f

（1）,f

（2）,jfc）;

构造函数子程序：

functiondelta=delta1（f,x,y）

a=f

（1）;

b=f

（2）;

c=f（3）;

delta=0;

fork=1:

15

delta=delta+（（x（k）+c）*y（k）-（a*x（k）+b））^2;

end

计算结果显示：

拟合出的方程为:

（x+-0.7110）y=11.2657x+-15.5024

均方差为:

0.33165089

总结：

指标选择，因题设方程

为非线性的，要转化为线性方程故需提指标为：

其驻点方程为:

计算结果显示：

4．设

，

，

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

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

（1）JACOBI迭代；

（1）Jacobi迭代

Jacobi迭代子程序：

function[x,k]=Jacobifl2（A,b,x0,eps,max1）

%jacobi按矩阵的分量算法

%x0初值，eps误差，max1最大迭代次数

n=length（x0）;

fork=1:

max1

x

（1）=（b

（1）-A（1,2:

n）*x0（2:

n））/A（1,1）;

fori=2:

n

x（i）=（b（i）-A（i,1:

i-1）*x0（1:

i-1）-A（i,i+1:

n）*x0（i+1:

n））/A（i,i）;

end

ifnorm（x'-x0,2）

return

else

x0=x';

end

end

Jacobi迭代计算程序：

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]';

x0=[000000]';

eps=10^（-4）;

max1=200;

[X,k]=Jacobifl2（A,b,x0,eps,max1）;

fprintf（'近似解x=\n[%4.8f\n%4.8f\n%4.8f\n%4.8f\n%4.8f\n%4.8f]\n迭代次数n=%d\n',X

（1）,X

（2）,X（3）,X（4）,X（5）,X（6）,k）;

计算结果显示：

近似解x=

[0.99998177

1.99995018

0.99997509

1.99995018

0.99997509

1.99996353]

迭代次数n=28

（2）GAUSS-SEIDEL迭代；

（2）GAUSS-SEIDEL迭代

GAUSS-SEIDEL迭代子程序：

function[x,k]=GSDD（A,b,x0,eps,max1）

n=length（x0）;

fork=1:

max1

x

（1）=（b

（1）-A（1,2:

n）*x0（2:

n）'）/A（1,1）;

fori=2:

n

x（i）=（b（i）-A（i,1:

i-1）*x（1:

i-1）'-A（i,i+1:

n）*x0（i+1:

n）'）/A（i,i）;

end

ifnorm（x-x0,2）

return

else

x0=x;

end

end

GAUSS-SEIDEL迭代计算程序：

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]';

x0=[000000];

eps=10^（-4）;

max1=200;

[X,k]=GSDD（A,b,x0,eps,max1）;

fprintf（'近似解x=\n[%4.8f\n%4.8f\n%4.8f\n%4.8f\n%4.8f\n%4.8f]\n迭代次数n=%d\n',X

（1）,X

（2）,X（3）,X（4）,X（5）,X（6）,k）;

计算结果显示：

近似解x=

[0.99996014

1.99995676

0.99996351

1.99996662

0.99997253

1.99998401]

迭代次数n=15

（3）SOR迭代（取

，找到迭代步数最少的

）。

（3）SOR迭代

SOR迭代子程序：

function[x,k,m]=SOR1（A,b,x0,eps,max1）

n=length（x0）;

a=x0;

fors=1:

19

x0=a;

w（s）=s/10;

x

（1）=x0

（1）+w（s）*（b

（1）-A（1,1:

n）*x0（1:

n）'）/A（1,1）;

fori=2:

n

x（i）=x0（i）+w（s）*（b（i）-A（i,1:

i-1）*x（1:

i-1）'-A（i,i:

n）*x0（i:

n）'）/A（i,i）;

end

k（s）=0;

whilenorm（x-x0,2）>=eps

x0=x;

x

（1）=x0

（1）+w（s）*（b

（1）-A（1,1:

n）*x0（1:

n）'）/A（1,1）;

fori=2:

n

x（i）=x0（i）+w（s）*（b（i）-A（i,1:

i-1）*x（1:

i-1）'-A（i,i:

n）*x0（i:

n）'）/A（i,i）;

end

k（s）=k（s）+1;

if（k（s）>=max1）

disp（'超出迭代次数上限！

\n'）;

return;

end

end

Data{s}=x;

end

[min1,index]=min（k）;

x=Data{index};

k=min1;

m=w（index）;

fprintf（'\n解x=\n%4.8f\n%4.8f\n%4.8f\n%4.8f\n%4.8f\n%4.8f',x）;

fprintf（'\n步数最小omega=%2.2f,步数n=%d。

\n',m,k）;

SOR迭代计算程序：

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]';

x0=[000000];

eps=10^（-4）;

max1=1000;

[X,k,m]=SOR1（A,b,x0,eps,max1）;

fprintf（'\n近似解x=\n%4.8f\n%4.8f\n%4.8f\n%4.8f\n%4.8f\n%4.8f',x）;

fprintf（'\n步数最小omega=%2.1f,步数n=%d\n',m,k）;

计算结果显示：

近似解x=

1.00000977

1.99999782

0.99999567

2.00001071

0.99999662

1.99999923

步数最小omega=1.20,步数n=9。

总结：

Jacobi迭代次数为28次；GAUSS-SEIDEL迭代次数为15次；SOR迭代次数在w=1.2时达到最小次数9次。

系数矩阵严格对角占优，则Jacobi迭代法与Gauss-Seidel迭代法均收敛；又因系数矩阵对称正定，且0＜w＜2，则SOR迭代法也收敛。

三种迭代法的结果分析及比较：

（1）Jacobi法，其设计思想是将线性方程组的求解归结为对角方程组求解过程的重复，计算规则简单而易于编写计算程序。

但收敛速度慢，谱半径相对偏大

（2）Gauss-Seidel充分利用新信息进行计算，其迭代的效果比Jacobi迭代好，谱半径比Jacobi法小。

（3）SOR迭代法，计算公式简单、编制程序容易，是求解大型稀疏方程组的一种有效方法，且若松弛因子选择得当超松弛法收敛速度比Gauss-Seidel迭代和Jacobi迭代都要快。

5．用逆幂迭代法求

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

。

逆幂迭代法子程序：

function[lamda,V,k]=NMSF2（A,v,eps,p）

[n,n]=size（A）;

A=A-p*eye（n）;

v0=v;

[tmax,tindex]=max（abs（v0））;

lamd0=v0（tindex）;

u0=v0/lamd0;

flag=0;

k=0;

while（flag==0）

V=A\u0;

[tmax,tindex]=max（abs（V））;

lamd1=V（tindex）;

u0=V/lamd1;

if（abs（（lamd0）^（-1）-（lamd1）^（-1）））<=eps

flag=1;

end

lamd0=lamd1;

k=k+1;

end（接下页）

lamda=（lamd1）^（-1）+p;

V=u0';

逆幂迭代计算程序：

A=[631;321;111];

v=[111]';

eps=0.001;

p=11;

[lamda,U,k]=NMSF2（A,v,eps,p）;

fprintf（'最接近11的特征值为:

%4.8f\n特征向量：

\n%4.8f\n%4.8f\n%4.8f\n',lamda,U

（1）,U

（2）,U（3））;

计算结果显示：

最接近于11的特征值为：

7.87313712

特征向量：

1.00000000

0.54922698

0.22545618

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

（1）

，

，

；

（2）

，

，

。

和精确解

比较，进行误差分析得到结论，图形显示精确解和数值解。

经典四阶R-K方法

R-K程序：

f=@（x,y）[-2*y

（1）+y

（2）+2*sin（x）;y

（1）-2*y

（2）+2*cos（x）-2*sin（x）];

g=@（x,y）[-2*y

（1）+y

（2）+2*sin（x）;998*y

（1）-999*y

（2）+

999*cos（x）-999*sin（x）];（接下页）

x=0:

0.1:

10;

Y1=2*exp（-x）+sin（x）;

Y2=2*exp（-x）+cos（x）;

x1=0:

10/10000:

10;

Y3=2*exp（-x1）+sin（x1）;

Y4=2*exp（-x1）+cos（x1）;

N=100;N1=10000;

[a1,b1]=ode23（f,[0:

0.1:

10],[23]）;

[a2,b2]=ode23（g,[0:

1/1000:

10],[23]）;

[a3,b3]=ode45（f,[0:

0.1:

10],[23]）;

[a4,b4]=ode45（g,[0:

1/1000:

10],[23]）;

n=length（b1）;

figure1=figure（'Color',[111]）;

dbt=plot（x,b1（:

1）,'r-',x,b1（:

2）,'k-',x,Y1,'g--',x,Y2,'y--'）;

legend（'y1','y2','精确解Y1','精确解Y2'）;

title（'四阶RK算法下的方程组的数值解&方程的精确解','FontSize',14,'FontWeight','Bold'）;

ylabel（'Y','FontSize',14,'FontWeight','Bold'）;

xlabel（'X','FontSize',14,'FontWeight','Bold'）;

set（dbt,'LineWidth',1.5）;gridon;

Q=b1（:

1）-Y1';W=b1（:

2）-Y2';

Q1=b2（:

1）-Y3';W1=b2（:

2）-Y4';

Q2=b3（:

1）-Y1';W2=b3（:

2）-Y2';

Q3=b4（:

1）-Y3';W3=b4（:

2）-Y4';

ER1=max（abs（Q））;ER2=max（abs（W））;

ER3=max（abs（Q1））;ER4=max（abs（W1））;

ER5=max（abs（Q2））;ER6=max（abs（W2））;

ER7=max（abs（Q3））;ER8=max（abs（W3））;

fprintf（'\n

（1）中y1在ode23下求出的最大误差为：

%4.8f\n

（1）中y2在ode23下求出的最大误差为：

%4.8f\n',ER1,ER2）;

fprintf（'\n

（2）中y1在ode23下求出的最大误差为：

%4.8f\n

（2）中y2在ode23下求出的最大误差为：

%4.8f\n',ER3,ER4）;

fprintf（'\n

（1）中y1在ode45下求出的最大误差为：

%4.8f\n

（1）中y2在ode45下求出的最大误差为：

%4.8f\n',ER5,ER6）;

fprintf（'\n

（2）中y1在ode45下求出的最大误差为：

%4.8f\n

（2）中y2在ode45下求出的最大误差为：

%4.8f\n',ER7,ER8）;

计算与显示程序：

（1）中y1在ode23下求出的最大误差为：

0.00169850

（1）中y2在ode23下求出的最大误差为：

0.00207422

（2）中y1在ode23下求出的最大误差为：

0.00000583

（2）中y2在ode23下求出的最大误差为：

0.00581329

（1）中y1在ode45下求出的最大误差为：

0.00052155

（1）中y2在ode45下求出的最大误差为：

0.00045793

（2）中y1在ode45下求出的最大误差为：

0.00000276

（2）中y2在ode45下求出的最大误差为：

0.00275331

总结：

对比2种方程在ode23，ode45命令下，与精确值的误差可以发现，在此题中，ode45的误差相对于ode23的误差小，较适合用在此类钢性问题中。

计算结果显示：

（1）Ode23的题1

（2）Ode3的题2

（3）

Ode45的题1

（4）Ode45的题2

7．用有限差分法求解边值问题（h=0.1），并