东南大学数值分析上机题作业MATLAB版.docx

东南大学数值分析上机题作业MATLAB版.docx

《东南大学数值分析上机题作业MATLAB版.docx》由会员分享，可在线阅读，更多相关《东南大学数值分析上机题作业MATLAB版.docx（29页珍藏版）》请在冰豆网上搜索。

东南大学数值分析上机题作业MATLAB版

东南大学-数值分析上机题作业-MATLAB版

D

>>P20T17

请输入N值:

10^2

精确值为：

0.740049

从大到小的顺序累加得SN=0.740049

从小到大的顺序累加得SN=0.740050

============================================================

>>P20T17

请输入N值:

10^4

精确值为：

0.749900

从大到小的顺序累加得SN=0.749852

1.3运行结果

从小到大的顺序累加得SN=0.749900

============================================================

>>P20T17

请输入N值:

10^6

精确值为：

0.749999

从大到小的顺序累加得SN=0.749852

从小到大的顺序累加得SN=0.749999

============================================================

1.4结果分析

按从大到小的顺序，有效位数分别为：

6，4，3。

按从小到大的顺序，有效位数分别为：

5，6，6。

可以看出，不同的算法造成的误差限是不同的，好的算法可以让结果更加精确。

当采用从大到小的顺序累加的算法时，误差限随着N的增大而增大，可见在累加的过程中，误差在放大，造成结果的误差较大。

因此，采取从小到大的顺序累加得到的结果更加精确。

2.Chapter2

2.1题目

（1）给定初值

及容许误差

，编制牛顿法解方程f（x）=0的通用程序。

（2）给定方程

易知其有三个根

由牛顿方法的局部收敛性可知存在

当

时，Newton迭代序列收敛于根x2*。

试确定尽可能大的

。

试取若干初始值，观察当

时Newton序列的收敛性以及收敛于哪一个根。

（3）通过本上机题，你明白了什么？

2.2程序

函数m文件：

fu.m

functionFu=fu（x）

Fu=x^3/3-x;

end

函数m文件：

dfu.m

functionFu=dfu（x）

Fu=x^2-1;

end

用Newton法求根的通用程序Newton.m

clear;

x0=input（'请输入初值x0:

'）;

ep=input（'请输入容许误差:

'）;

flag=1;

whileflag==1

x1=x0-fu（x0）/dfu（x0）;

ifabs（x1-x0）

flag=0;

end

x0=x1;

end

fprintf（'方程的一个近似解为：

%f\n',x0）;

寻找最大δ值的程序：

Find.m

clear

eps=input（'请输入搜索精度：

'）;

ep=input（'请输入容许误差：

'）;

flag=1;

k=0;

x0=0;

whileflag==1

sigma=k*eps;

x0=sigma;

k=k+1;

m=0;

flag1=1;

whileflag1==1&&m<=10^3

x1=x0-fu（x0）/dfu（x0）;

ifabs（x1-x0）

flag1=0;

end

m=m+1;

x0=x1;

end

ifflag1==1||abs（x0）>=ep

flag=0;

end

end

fprintf（'最大的sigma值为：

%f\n',sigma）;

2.3运行结果

（1）寻找最大的

值。

算法为：

将初值x0在从0开始不断累加搜索精度eps，带入Newton迭代公式，直到求得的根不再收敛于0为止，此时的x0值即为最大的sigma值。

运行Find.m，得到在不同的搜索精度下的最大sigma值。

>>Find

请输入搜索精度：

10^-6

请输入容许误差：

10^-6

最大的sigma值为：

0.774597

>>Find

请输入搜索精度：

10^-4

请输入容许误差：

10^-6

最大的sigma值为：

0.774600

>>Find

请输入搜索精度：

10^-2

请输入容许误差：

10^-6

最大的sigma值为：

0.780000

（2）运行Newton.m

在

内取初值，运行结果如下：

X0

Xk

-1000

-1.732051

-500

-1.732051

-100

-1.732051

-10

-1.732051

-5

-1.732051

-2.5

-1.732051

-1.5

-1.732051

可见，在

区间内取初值，Newton序列收敛，且收敛于根

。

在

内取初值，运行结果如下：

X0

Xk

-0.95

1.732051

-0.85

1.732051

-0.8

1.732051

-0.774598

1.732051

可见，在

内取初值，Newton序列收敛，且收敛于根

。

在

内内取初值，运行结果如下：

X0

Xk

-0.774596

0.000000

-0.55

0.000000

-0.35

0.000000

-0.15

0.000000

0.05

0.000000

0.25

0.000000

0.45

0.000000

0.65

0.000000

0.774596

0.000000

可见，在

内取初值，Newton序列收敛，且收敛于根0。

在

内取初值，运行结果如下：

X0

Xk

0.774598

-1.732051

0.8

-1.732051

0.85

-1.732051

0.95

-1.732051

可见，在

内取初值，Newton序列收敛，且收敛于根

在

内取初值，运行结果如下：

X0

Xk

1.5

1.732051

2.5

1.732051

5

1.732051

10

1.732051

100

1.732051

500

1.732051

1000

1.732051

可见，在

内取初值，Newton序列收敛，且收敛于根

3.Chapter3

3.1题目

对于某电路的分析，归结为求解线性方程组RI=V，其中

（1）编制解n阶线性方程组

的列主元高斯消去法的通用程序；

（2）用所编程序线性方程组

，并打印出解向量，保留5位有效数字；

（3）本题编程之中，你提高了哪些编程能力？

3.2程序

n=input（'请输入线性方程组阶数:

n='）;

b=zeros（1,n）;

A=input（'请输入系数矩阵：

A=\n'）;

b（1,:

）=input（'请输入线性方程组右端向量：

b=\n'）;

b=b';

C=[A,b];

fori=1:

n-1

[maximum,index]=max（abs（C（i:

n,i）））;

index=index+i-1;

T=C（index,:

）;

C（index,:

）=C（i,:

）;

C（i,:

）=T;

fork=i+1:

n

ifC（k,i）~=0

C（k,:

）=C（k,:

）-C（k,i）/C（i,i）*C（i,:

）;

end

end

end

%%回代求解

x=zeros（n,1）;

x（n）=C（n,n+1）/C（n,n）;

fori=n-1:

-1:

1

x（i）=（C（i,n+1）-C（i,i+1:

n）*x（i+1:

n,1））/C（i,i）;

end

disp（'方程组的解为:

'）;

fprintf（'%.5g\n',x）;

3.3运行结果

运行程序，输入系数矩阵和方程组右端列向量。

运行过程与结果如下图所示：

>>P126T39

请输入线性方程组阶数:

n=4

请输入系数矩阵：

A=

[136.0190.8600;90.8698.81-67.590;0-67.59132.0146.26;0046.26177.17]

请输入线性方程组右端向量：

b=

[-33.25449.7928.067-7.324]

方程组的解为:

-2957.4

4426.6

2495

-651.49

>>P126T39

请输入线性方程组阶数:

n=9

请输入系数矩阵：

A=

[31-13000-10000;-1335-90-110000;0-931-1000000;00-1079-30000-9;000-3057-70-50;0000-747-3000;00000-304100;0000-50027-2;000-9000-229]

请输入线性方程组右端向量：

b=

[-1527-230-2012-7710]

方程组的解为:

-0.28923

0.34544

-0.71281

-0.22061

-0.4304

0.15431

-0.057823

0.20105

0.29023

可看出，算得的该线性方程组的解向量为：

[-0.289230.34544-0.71281-0.22061-0.43040.15431-0.0578230.201050.29023]

4.Chapter4

4.1题目

（1）编制求第一型3次样条插值函数的通用程序；

（2）已知汽车门曲线型值点的数据如下：

i

0

1

2

3

4

5

6

7

8

9

10

Xi

0

1

2

3

4

5

6

7

8

9

10

Yi

2.51

3.30

4.04

4.70

5.22

5.54

5.78

5.40

5.57

5.70

5.80

端点条件为

，用所编程序求车门的3次样条插值函数S（x）,并打印出S（i+0.5），i=0,1，…，9。

4.2程序

clear

digits（6）;

n=input（'请输入节点数：

n='）;

xn=zeros（1,n）;

yn=zeros（1,n）;

xn（1,:

）=input（'请输入节点坐标：

'）;

yn（1,:

）=input（'请输入节点处函数值：

'）;

dy0=input（'请输入左边界条件：

y’（x0）='）;

dyn=input（'请输入右边界条件：

y’（xn）='）;

%====================求d====================%

d=zeros（n,1）;

h=zeros（1,n-1）;

f1=zeros（1,n-1）;

f2=zeros（1,n-2）;

fori=1:

n-1

h（i）=xn（i+1）-xn（i）;

f1（i）=（yn（i+1）-yn（i））/h（i）;

end

fori=2:

n-1

f2（i）=（f1（i）-f1（i-1））/（xn（i+1）-xn（i-1））;

d（i）=6*f2（i）;

end

d（i）=6*（f1

（1）-dy0）/h

（1）;

d（n）=6*（dyn-f1（n-1））/h（n-1）;

%====================求Mi====================%

A=zeros（n）;

miu=zeros（1,n-2）;

lamda=zeros（1,n-2）;

fori=1:

n-2

miu（i）=h（i）/（h（i）+h（i+1））;

lamda（i）=1-miu（i）;

end

A（1,2）=1;

A（n,n-1）=1;

fori=1:

n

A（i,i）=2;

end

fori=2:

n-1

A（i,i-1）=miu（i-1）;

A（i,i+1）=lamda（i-1）;

end

M=A\d;

%====================回代求插值函数====================%

symsx;

fori=1:

n-1;

Sx（i）=collect（yn（i）+（f1（i）-（M（i）/3+M（i+1）/6）*h（i））*（x-xn（i））+M（i）/2*（x-xn（i））^2+（M（i+1）-M（i））/（6*h（i））*（x-xn（i））^3）;

Sx（i）=vpa（Sx（i）,6）;

end

S=zeros（1,n-1）;

fori=1:

n-1

x=xn（i）+0.5;

S（i）=yn（i）+（f1（i）-（M（i）/3+M（i+1）/6）*h（i））*（x-xn（i））+M（i）/2*（x-xn（i））^2+（M（i+1）-M（i））/（6*h（i））*（x-xn（i））^3;

end

%====================打印结果====================%

disp（'S（x）='）;

fori=1:

n-1

formatshort;

fprintf（'%s（%d

disp（'======================================================================'）;

end

disp（'S（i+0.5）'）

disp（'ix（i+0.5）S（i+0.5）'）;

fori=1:

n-1

fprintf（'%d%.5f%.5f\n',i,xn（i）+0.5,S（i））;

end

4.3运行结果

>>P195T37

请输入节点数：

n=11

请输入节点坐标：

[012345678910]

请输入节点处函数值：

[2.513.304.044.705.225.545.785.405.575.705.80]

请输入左边界条件：

y’（x0）=0.8

请输入右边界条件：

y’（xn）=0.2

S（x）=

0.79*x+0.0158344*x^2-0.0158344*x^3+2.51（0

======================================================================

0.830013*x-0.0241785*x^2-0.00249676*x^3+2.49666（1

======================================================================

0.809832*x-0.0140879*x^2-0.00417854*x^3+2.51012（2

======================================================================

0.315407*x^2-0.178653*x-0.0407891*x^3+3.4986（3

======================================================================

6.9313*x-1.46208*x^2+0.107335*x^3-5.98133（4

======================================================================

4.1762*x^2-21.2601*x-0.26855*x^3+41.0043（5

======================================================================

53.8449*x-8.3413*x^2+0.426866*x^3-109.206（6

======================================================================

6.27011*x^2-48.435*x-0.268915*x^3+129.447（7

======================================================================

14.4854*x-1.59494*x^2+0.0587951*x^3-38.3403（8

======================================================================

13.2458*x-1.45831*x^2+0.053735*x^3-34.5615（9

======================================================================

S（i+0.5）

ix（i+0.5）S（i+0.5）

10.500002.90698

21.500003.67885

32.500004.38136

43.500004.98822

54.500005.38326

65.500005.72372

76.500005.59435

87.500005.43012

98.500005.65892

109.500005.73172

5.Chapter5

5.1题目

用Romberg求积法计算积分

的近似值，要求误差不超过

。

5.2程序

%被积函数m文件：

fx.m

functionFx=fx（x）

Fx=1/（1+100*x*x）;

end

%Romberg求积法计算积分的通用程序

functionRomberg（）

clear;

a=input（'请输入积分下限：

a='）;

b=input（'请输入积分上限：

b='）;

eps=input（'请输入允许精度：

eps='）;

%========计算Tn========%

functionTn=T（n）

Tn=0;

h=（b-a）/n;

x=zeros（1,n+1）;

fork=1:

n+1

x（k）=a+（k-1）*h;

end

forj=1:

n

Tn=Tn+h*（fx（x（j））+fx（x（j+1）））/2;

end

end

%========计算Sn========%

functionSn=S（n）

Sn=4/3*T（2*n）-1/3*T（n）;

end

%========计算Cn========%

functionCn=C（n）

Cn=16/15*S（2*n）-1/15*S（n）;

end

%========计算Rn========%

functionRn=R（n）

Rn=64/63*C（2*n）-1/63*C（n）;

end

%========计算满足允许精度的Rn，并打印输出========%

i=1;

flag=1;

whileflag==1

ifabs（R（2^i）-R（2^（i-1）））/255

flag=0;

end

i=i+1;

end

fprintf（'该积分的值为：

%f\n',R（2^（i-1）））;

end

5.3运行结果

>>Romberg

请输入积分下限：

a=-1

请输入积分上限：

b=1

请输入允许精度：

eps=0.5*10^-7

该积分的值为：

0.2942255

5.4结果分析

手动化简该定积分并最终求得的值为：

0.294225534860747，误差限为：

，可见，程序完成了计算要求。

6.Chapter6

6.1题目

常微分方程初值问题数值解

（1）编制RK4方法的通用程序；

（2）编制AB4方法的通用程序（由RK4提供初值）；

（3）编制AB4-AM4预测校正方法通用程序（由RK4提供初值）；

（4）编制带改进的AB4-AM4预测校正方法通用程序（由RK4提供初值）；

（5）对于初值问题

取步长h=0.1，应用

（1）-（4）中的四种方法进行计算，并将计算结果和精确解

作比较；

（6）通过本上机题，你能得到哪些结论？

6.2程序

%f（x，y）函数m文件：

fxy.m

functionFXY=fxy（x,y）

FXY=-x*x*y*y;

end

%精确解y（x）函数m文件：

fx.m

functionFX=fx（x）

FX=3/（1+x*x*x）;

end

%RK4法通用程序

functionRK4（）

clear;

x

（1）=input（'请输入初始x值：

x0='）;

y

（1）=input（'请输入初值条件：

y（x0）='）;

N=input（'请输入计算步长：

N='）;

h=input（'请输入步长：

h='）;

fori=1:

N-1

x（i+1）=x（i）+h;

k1=fxy（x（i）,y（i））;

k2=fxy（x（i）+0.5*h,y（i）+0.5*h*k1）;

k3=fxy（x（i）+0.5*h,y（i）+0.5*h*k2）;

k4=fxy（x（i）+h,y（i）+h*k3）;

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

end

disp（'ixiyiy（xi）y（xi）-yi'）;

disp（'-------------------------------------------------------------'）;

fori=1:

N

fprintf（'%d%f%f%f%f\n',i,x（i）,y（i）,fx（x（i））,fx（x（i））-y（i））;

disp（'-------------------------------------------------------------'）;

end

end

%AB4法通用程序

functionAB4（）

clear;

x

（1）=input（'请输入初始x值：

x0='）;

y

（1）=input（'请输入初值条件：

y（x0）='）;

N=input（'请输入计算步长：

N='）;

h=input（'请输入步长：

h='）;

fori=1:

N-1

x（i+1）=x（i）+h;

k1=fxy（x（i）,y（i））;

k2=fxy（x（i）+0.5*h,y（i）+0.5*h*k1）;

k3=fxy（x（i）+0.5*h,y（i）+0.5*h*k2）;

k4=fxy（x（i）+h,y（i）+h*k3）;

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

end

fori=4:

N-1

y（i+1）=y（i）+h/24*（55*fxy（x（i）,y（i））-59*fxy（x（i-1）,y（i-1））+37*fx