收藏
下载资源下载资源 加入VIP,免费下载

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

上传人：b****8 文档编号：10738576 上传时间：2023-02-22 格式：DOCX 页数：27 大小：155.67KB
下载 相关 举报
数值计算方法上机实习题NEW.docx_第1页
第1页 / 共27页
数值计算方法上机实习题NEW.docx_第2页
第2页 / 共27页
数值计算方法上机实习题NEW.docx_第3页
第3页 / 共27页
数值计算方法上机实习题NEW.docx_第4页
第4页 / 共27页
数值计算方法上机实习题NEW.docx_第5页
第5页 / 共27页
点击查看更多>>
下载资源
资源描述

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

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

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

数值计算方法上机实习题NEW

数值计算方法上机实习题

1．设

，

（1）由递推公式

，从

出发，计算

；

程序如下：

functionI=myhs（I0,n）

ifn>=1

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

elseifn==0

I=I0;

end

命令行窗口输入：

I0=0.1822;

I1=myhs（I0,20）;

I0=0.1823;

I2=myhs（I0,20）;

运行结果：

当I0=0.1822时，I20=-1.1593e+10。

当I0=0.1823时，I20=-2.0558e+9。

（2）

，

用

，计算

；

程序如下：

functionI=myhs2（I20,n）

if（n<20）

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

elseifn==20

I=I20;

end

命令行窗口输入：

I20=0;

I1=myhs2（I20,20）;

I20=10000;

I2=myhs2（I20,20）;

运行结果：

当I20=0时，I0=0.182********3955。

当I20=10000时，I0=0.182********8812。

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

根据上式，

假设In的真值为I*，误差为en，即e=I*-In。

综合递推式，有en=-5*en-1，这意味着哪怕开始只有一点点误差，只要n足够大，按照这种计算一步误差增长5倍的方式，所得的结果总是不可信的，因此整个算法是不稳定的。

根据上式，假设In的真值为I*，误差为en，即e=I*-In。

综合递推式，有en-1.=（-1/5）*en，按照这种计算误差会以每步缩小到1/5的方式进行，根据

（2）得到的结果而言，该算法是相对稳定的。

2．求方程

的近似根，要求

，并比较计算量。

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

输入程序如下：

a=0;

b=1;

n=0;

while（abs（b-a）>5*1e-4）

c=（a+b）/2;

y=exp（c）+10*c-2;

if（y>0）

b=c;

else

a=c;

end

n=n+1;

end

c

得到结果：

c=0.090332，n=11

（2）取初值

，并用迭代

；

输入程序如下：

xk=0;

xki=（2-exp（xk））/10;

n=1;

while（abs（xki-xk）>5*1e-4）

xk=xki;

xki=（2-exp（xk））/10;

n=n+1;

end

xki

程序运行结果为：

xki=0.090513，n=4

（3）加速迭代的结果；

输入程序如下

xk=0;

yk=（2-exp（xk））/10;

zk=（2-exp（yk））/10;

xki=xk-（yk-xk）^2/（zk-2*yk+xk）;

n=1;

while（abs（xki-xk）>5*1e-4）

xk=xki;

yk=（2-exp（xk））/10;

zk=（2-exp（yk））/10;

xki=xk-（yk-xk）^2/（zk-2*yk+xk）;

n=n+1;

end

xki

程序运行的结果为

xki=0.090525，n=2

（4）取初值

，并用牛顿迭代法；

输入程序如下：

xk=0;

xki=xk-（exp（xk）+10*xk-2）/（exp（xk）+10）;

n=1;

while（abs（xki-xk）>5*1e-4）

xk=xki;

xki=xk-（exp（xk）+10*xk-2）/（exp（xk）+10）;

n=n+1;

end

xki

程序运行的结果为

xki=0.090525，n=2

（5）分析绝对误差。

通过程序x=solve（'exp（x）+10*x-2=0'）;

解得x=0.090525

用二分法求得x的值为0.090332，迭代次数为11，绝对误差为0.000193。

用不动点迭代法求得x的值为0.090513，迭代次数为4，绝对误差为0.000012。

用加速迭代法和牛顿迭代法求得x的值均为0.090525，迭代次数均为2，绝对误差为0。

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=（ax+b）/（c+x），将原模型变为（c+x）y=（ax+b），采用非线性最小二乘法。

应选取参数a、b、c使得

达到极小值。

具体方法：

对Z分别求a、b、c的偏导，并使偏导为0，相应的方程组如下

对上述方程组进行整理，写成AX=B的形式，利用MATLAB求解X，对应的就是a、b、c三个参数的值。

具体程序如下：

%方程组求解和画拟合曲线

C=xlsread（'test3.xlsx'）;

x=C（1,:

）;

y=C（2,:

）;

A=[sum（x.^2）,sum（x）,-sum（x.*y）;sum（x）,15,-sum（y）;sum（x.*y）,sum（y）,-sum（y.^2）];

B=[sum（x.^2.*y）;sum（x.*y）;sum（x.*y.^2）];

Z=A\B;

a=Z

（1）;

b=Z

（2）;

c=Z（3）;

X=linspace（1,50,50）;

Y=zeros（1,50）;

fori=1:

50

Y（i）=（a*X（i）+b）/（c+X（i））;

end

plot（x,y,'*',X,Y,'-'）;

%计算方差

z1=0;

fori=2:

16

z1=z1+Y（i）;

end

z1=z1/15;

sum=0;

fori=2:

16

sum=sum+（Y（i）-z1）^2;

end

fc=sqrt（sum/15）;

运行结果为：

a=11.3400，b=-12.4953，c=-0.3403

拟合方程为y=（11.34x-12.4953）/（-0.3403+x）,均方差为1.2285

4．设

，

，

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

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

（1）JACOBI迭代；

解：

输入程序如下

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

D=diag（diag（A））;

L=D-tril（A）;

U=D-triu（A）;

Bj=D\（L+U）;

[V,landa]=eig（Bj）;

Maxlanda=abs（landa（1,1））;

fori=2:

size（landa）

if（abs（landa（i,i））>Maxlanda）

Maxlanda=abs（landa（i,i））;

end

end

运行结果为：

Maxlanda=0.6830

即雅克比迭代矩阵半径=0.683<1，雅克比迭代收敛。

%自定义的雅克比迭代函数

function[x,iter]=jacobi（A,b,x0,tol）

D=diag（diag（A））;

L=D-tril（A）;

U=D-triu（A）;

x=D\（b+L*x0+U*x0）;

iter=1;

error=norm（x-x0）;

while（error>0.0001）

fprintf（'%d%f%f%f%f%f%f%f\n',iter,x

（1）,x

（2）,x（3）,x（4）,x（5）,x（6）,error）;

if（error

break;

end

x0=x;

x=D\（b+L*x0+U*x0）;

error=norm（x-x0）;

iter=iter+1;

end

fprintf（'%d%f%f%f%f%f%f%f\n',iter,x

（1）,x

（2）,x（3）,x（4）,x（5）,x（6）,error）;

命令行窗口输入以下程序即可得到近似解：

formatshort;

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

tol=0.0001;

x0=zeros（size（B））;

jacobi（A,B,x0,tol）;

运行结果为：

10.0000001.250000-0.5000001.250000-0.5000001.5000002.423840

20.6250001.0000000.5000001.0000000.5000001.2500001.605654

30.5000001.6562500.3125001.6562500.3125001.7500001.094196

40.8281251.5312500.7656251.5312500.7656251.6562500.747228

50.7656251.8398440.6796881.8398440.6796881.8828130.510360

60.9199221.7812500.8906251.7812500.8906251.8398440.348582

70.8906251.9252930.8505861.9252930.8505861.9453130.238086

80.9626461.8979490.9489751.8979490.9489751.9252930.162616

90.9489751.9651490.9302981.9651490.9302981.9744870.111069

100.9825741.9523930.9761961.9523930.9761961.9651490.075861

110.9761961.9837420.9674841.9837420.9674841.9880980.051814

120.9918711.9777910.9888951.9777910.9888951.9837420.035390

130.9888951.9924150.9848311.9924150.9848311.9944480.024172

140.9962081.9896390.9948201.9896390.9948201.9924150.016510

150.9948201.9964620.9929231.9964620.9929231.9974100.011276

160.9982311.9951670.9975831.9951670.9975831.9964620.007702

170.9975831.9983490.9966991.9983490.9966991.9987920.005260

180.9991751.9977450.9988731.9977450.9988731.9983490.003593

190.9988731.9992300.9984601.9992300.9984601.9994360.002454

200.9996151.9989480.9994741.9989480.9994741.9992300.001676

210.9994741.9996410.9992821.9996410.9992821.9997370.001145

220.9998201.9995090.9997551.9995090.9997551.9996410.000782

230.9997551.9998320.9996651.9998320.9996651.9998770.000534

240.9999161.9997710.9998861.9997710.9998861.9998320.000365

250.9998861.9999220.9998441.9999220.9998441.9999430.000249

260.9999611.9998930.9999471.9998930.9999471.9999220.000170

270.9999471.9999640.9999271.9999640.9999271.9999730.000116

280.9999821.9999500.9999751.9999500.9999751.9999640.000079

线性方程组的雅克比迭代近似解为：

x=[0.9999821.9999500.9999751.9999500.9999751.999964]

误差为0.000079

（2）GAUSS-SEIDEL迭代；

解：

首先，求GAUSS-SEIDEL迭代的收敛性，输入如下程序：

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

D=diag（diag（A））;

L=D-tril（A）;

U=D-triu（A）;

Bg=（D+L）\U;

[V,landa]=eig（Bg）;

Maxlanda=abs（landa（1,1））;

fori=2:

size（landa）

if（abs（landa（i,i））>Maxlanda）

Maxlanda=abs（landa（i,i））;

end

end

运行结果为

Maxlanda=0.4806

即GAUSS-SEIDEL迭代矩阵半径=0.4806<1，GAUSS-SEIDEL迭代收敛。

接下来，求GAUSS-SEIDEL迭代的近似解：

%自定义的GAUSS-SEIDEL迭代函数

function[x,iter]=gauseidel（A,b,x0,tol）

D=diag（diag（A））;

L=D-tril（A）;

U=D-triu（A）;

x=（D+L）\（b+U*x0）;

iter=1;

error=norm（x-x0）;

while（error>0.0001）

fprintf（'%d%f%f%f%f%f%f%f\n',iter,x

（1）,x

（2）,x（3）,x（4）,x（5）,x（6）,error）;

if（error

break;

end

x0=x;

x=（D+L）\（b+U*x0）;

error=norm（x-x0）;

iter=iter+1;

end

fprintf（'%d%f%f%f%f%f%f%f\n',iter,x

（1）,x

（2）,x（3）,x（4）,x（5）,x（6）,error）;

命令行窗口输入如下程序：

formatshort;

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

tol=0.0001;

x0=zeros（size（B））;

gauseidel（A,B,x0,tol）;

运行结果为：

10.0000001.250000-0.8125001.453125-1.1757811.9970703.115282

20.6757810.5839840.2165530.732971-0.3299711.5283551.847408

30.3292391.139336-0.2195021.140073-0.6877641.7268160.975595

40.5698520.880720-0.0034580.936461-0.5225911.6315120.499274

50.4542951.004914-0.1092351.033087-0.6016221.6777140.240174

60.5095000.944911-0.0585270.986851-0.5635121.6555100.115337

70.4829400.973755-0.0828491.009099-0.5818361.6661710.055438

80.4957140.959900-0.0711580.998402-0.5730331.6610480.026645

90.4895760.966559-0.0767771.003542-0.5772631.6635100.012805

100.4925250.963359-0.0740771.001072-0.5752301.6623270.006154

110.4911080.964896-0.0753741.002259-0.5762071.6628950.002957

120.4917890.964157-0.0747511.001689-0.5757381.6626220.001421

130.4914620.964513-0.0750501.001963-0.5759631.6627530.000683

140.4916190.964342-0.0749061.001831-0.5758551.6626900.000328

150.4915430.964424-0.0749761.001894-0.5759071.6627210.000158

160.4915800.964384-0.0749421.001864-0.5758821.6627060.000076

线性方程组的雅克比迭代近似解为：

x=[0.4915800.964384-0.0749421.001864-0.5758821.662706]

误差为0.000076

（3）SOR迭代（取

，找到迭代步数最少的

）。

解：

首先看SOR迭代的收敛性，输入以下程序：

formatshort;

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

D=diag（diag（A））;

L=D-tril（A）;

U=D-triu（A）;

w=0.1:

0.1:

1.9;

n=length（w）;

l=1;

fori=1:

n

Bs=（D+w（i）*L）\（D-w（i）*D+w（i）*U）;

[V,landa]=eig（Bs）;

Maxlanda=abs（landa（1,1））;

forj=2:

size（landa）

if（abs（landa（j,j））>Maxlanda）

Maxlanda=abs（landa（j,j））;

end

end

if（Maxlanda<1）

W（l）=w（i）;

l=l+1;

end

end

fori=length（W）

fprintf（'%f\n',W（i））;

end

运行结果为：

W=[0.10000.20000.30000.40000.50000.60000.70000.80000.90001.00001.1000]

即SOR迭代的w*为W内的值时，该迭代方法是收敛的。

然后求SOR迭代次数最少时的w*，输入以下程序：

%自定义的sor迭代函数

function[x,iter,error]=sor（A,b,x0,w,tol）

D=diag（diag（A））;

L=D-tril（A）;

U=D-triu（A）;

x=（D-w*L）\（D*x0+w*（b-D*x0+U*x0））;

iter=1;

error=norm（x-x0）;

while（error>tol）

x0=x;

x=（D-w*L）\（D*x0+w*（b-D*x0+U*x0））;

error=norm（x-x0）;

iter=iter+1;

end

在命令行输入程序：

formatshort;

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

tol=0.0001;

x0=zeros（size（B））;

w=[0.10.20.30.40.50.60.70.80.91.01.1];

[X,Iter,Error]=sor（A,B,x0,0.1,tol）;

fori=2:

length（w）

x0=zeros（size（B））;

[x,iter,error]=sor（A,B,x0,w（i）,tol）;

if（Iter>iter）

X=x;

Iter=iter;

W=w（i）;

Error=error;

end

end

运行得到：

Iter=12X=[1.00002.00001.00002.00001.00002.0000]

Error=6.8688e-05W=1.1000

即SOR迭代收敛的条件下，当w*=1.1时，迭代次数最少为12，对应的近似解为x=[1.00002.00001.00002.00001.00002.0000]，误差为6.8688e-05

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

（1）

，

，

；

（2）

，

，

。

和精确解

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

解：

1）利用欧拉公式求解式

（1）的数值解，然后将得到的数值解与精确解通过图形显示，步长设置为0.1，输入如下程序：

clear;

a=0;b=10;N=100;y1（1,:

）=[2,3];

h=（b-a）/N;x=a:

h:

b;

forn=1:

N

y1（n+1,1）=y1（n,1）+h*（-2*y1（n,1）+y1（n,2

展开阅读全文
相关资源
猜你喜欢
相关搜索

当前位置：首页 > 初中教育 > 初中作文

copyright@ 2008-2022 冰豆网网站版权所有

经营许可证编号:鄂ICP备2022015515号-1