东南大学《数值分析报告》上机题.docx

东南大学《数值分析报告》上机题.docx

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

东南大学《数值分析报告》上机题

数值分析上机题1

设

，其精确值为

。

（1）编制按从大到小的顺序

，计算

的通用程序。

（2）编制按从小到大的顺序

，计算

的通用程序。

（3）按两种顺序分别计算

，

，

，并指出有效位数。

（编制程序时用单精度）

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

程序代码（matlab编程）：

clc

clear

a=single（1./（[2:

10^7].^2-1））;

S1

（1）=single（0）;

S1

（2）=1/（2^2-1）;

forN=3:

10^2

S1（N）=a

（1）;

fori=2:

N-1

S1（N）=S1（N）+a（i）;

end

end

S2

（1）=single（0）;

S2

（2）=1/（2^2-1）;

forN=3:

10^2

S2（N）=a（N-1）;

fori=linspace（N-2,1,N-2）

S2（N）=S2（N）+a（i）;

end

end

S1表示按从大到小的顺序的SN

S2表示按从小到大的顺序的SN

计算结果

从大到小的顺序的值

从小到大的顺序的值

精确值

有效位数

从大到小

从小到大

0.740049

0.74005

0.740049

6

5

0.749852

0.7499

0.7499

4

4

0.749852

0.749999

0.749999

3

6

通过本上机题，看出按两种不同的顺序计算的结果是不相同的，按从大到小的顺序计算的值与精确值有较大的误差，而按从小到大的顺序计算的值与精确值吻合。

从大到小的顺序计算得到的结果的有效位数少。

计算机在进行数值计算时会出现“大数吃小数”的现象，导致计算结果的精度有所降低，我们在计算机中进行同号数的加法时，采用绝对值较小者先加的算法，其结果的相对误差较小。

数值分析上机题2

20．（上机题）Newton迭代法

（1）给定初值

及容许误差

，编制Newton法解方程

根的通用程序。

（2）给定方程

，易知其有三个根

，

，

。

1．由Newton方法的局部收敛性可知存在

，当

时，Newton迭代序列收敛于根

。

试确定尽可能大的

。

2．试取若干初始值，观察当

，

，

，

，

时Newton序列是否收敛以及收敛于哪一个根。

MATLAB程序

问题1

clc

clear

dx=0.5;

x

（1）=0.5;

while（dx>1e-6）

i=1;

error=1;

while（error>1e-8）

x（i+1）=x（i）-（1/3*x（i）^3-x（i））/（x（i）^2-1）;

error=abs（x（i+1）-x（i））;

i=i+1;

end

if（x（i）==0）

x

（1）=x

（1）+dx;

else

dx=dx/2;

x

（1）=x

（1）-dx;

end

end

经计算，最大的

为0.774596

问题2

clc

clear

x2

（1）=1e14;

i=1;

error=1;

while（error>1e-8）

x2（i+1）=x2（i）-（1/3*x2（i）^3-x2（i））/（x2（i）^2-1）;

error=abs（x2（i+1）-x2（i））;

i=i+1;

if（i>1e4）

break

end

end

对于不同得初始值收敛于不同的根,

在（-∞，-1）收敛于

，在（-0.774，0.774）收敛于0，在（1，+∞）收敛于

，但在（0.774，1）和（－1，0.774）均可能收敛于

和

。

分析：

对于不同的初值，迭代序列会收敛于不同的根，所以在某个区间求根对于初值的选取有很大的关系。

产生上述结果的原因是区间不满足大围收敛的条件。

数值分析上机题3

39．（上机题）列主元三角分解法对于某电路的分析，归结为求解线性方程组RI=V。

（1）编制解n阶线性方程组Ax=b的列主元三角分解法的通用程序；

（2）用所编制的程序解线性方程组RI=V，并打印出解向量，保留五位有效数；

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

程序：

clc

clear

A=[31,-13,0,0,0,-10,0,0,0

-13,35,-9,0,-11,0,0,0,0

0,-9,31,-10,0,0,0,0,0

0,0,-10,79,-30,0,0,0,-9

0,0,0,-30,57,-7,0,-5,0

0,0,0,0,-7,47,-30,0,0

0,0,0,0,0,-30,41,0,0

0,0,0,0,-5,0,0,27,-2

0,0,0,-9,0,0,0,-2,29];

b=[-15,27,-23,0,-20,12,-7,7,10]';

[m,n]=size（A）;

Ap=[A,b];

x=zeros（n,1）;

fori=1:

m-1

j=i;

[maxa,maxi]=max（abs（Ap（i:

end,j）））;

maxi=maxi+i-1;

if（maxa~=0）

mid=Ap（maxi,:

）;

Ap（maxi,:

）=Ap（i,:

）;

Ap（i,:

）=mid;

fork=i:

m

Ap（i+1:

m,:

）=Ap（i+1:

m,:

）-Ap（i+1:

m,j）*（Ap（i,:

）./maxa）;

end

end

end

fori=linspace（m,1,m）

x（i）=（Ap（i,end）-Ap（i,1:

end-1）*x）/Ap（i,i）;

end

结果：

方程的解为（保留5位有效数字）：

x1=-0.28923，x2=0.34544，x3=-0.71281，

x4=-0.22061，x5=-0.43040，x6=0.15431，

x7=-0.057823，x8=0.20105，x9=0.29023。

习题4

37.（上机题）3次样条插值函数

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

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

0

1

2

3

4

5

6

7

8

9

10

2.51

3.30

4.04

4.70

5.22

5.54

5.78

5.40

5.57

5.70

5.80

端点条件为

=0.8，

=0.2。

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

程序：

（1）

clc

clear

%%

x=[0,1,2,3,4,5,6,7,8,9,10];

y=[2.51,3.30,4.04,4.7,5.22,5.54,5.78,5.40,5.57,5.70,5.80];

y1=0.8;

yend=0.2;

%%___________________________________________

n=size（x,2）-1;

h=x（2:

end）-x（1:

end-1）;

miu=h（1:

end-1）./（h（1:

end-1）+h（2:

end））;

lamda=1-miu;

f1=[y1,（y（2:

end）-y（1:

end-1））./h,yend];%f[xn-1,xn]

f2=[f1（2:

end）-f1（1:

end-1）]./[h

（1）,h（1:

end-1）+h（2:

end）,h（end）];%f[xn-1,xn,xn+1]

A=2.*eye（n+1）;

A（2:

end,1:

end-1）=A（2:

end,1:

end-1）+diag（[miu,1]'）;

A（1:

end-1,2:

end）=A（1:

end-1,2:

end）+diag（[1,lamda]'）;

M=A\（6*f2'）;

Sx=[y（1:

end-1）',（（y（2:

end）-y（1:

end-1））./h）'-（（1/3*M（1:

end-1）+1/6*M（2:

end））.*h'）,1/2*M（1:

end-1）,1/6*（M（2:

end）-M（1:

end-1））./h'];

%%

xx=input（’x=’）;

forj=2:

n+1

ifxx

S=Sx（j-1,:

）*[1,xx-x（j-1）,（xx-x（j-1））^2,（xx-x（j-1））^3]';

break

end

end

（2）

clc

clear

%%

x=[0,1,2,3,4,5,6,7,8,9,10];

y=[2.51,3.30,4.04,4.7,5.22,5.54,5.78,5.40,5.57,5.70,5.80];

y1=0.8;

yend=0.2;

%%___________________________________________

n=size（x,2）-1;

h=x（2:

end）-x（1:

end-1）;

miu=h（1:

end-1）./（h（1:

end-1）+h（2:

end））;

lamda=1-miu;

f1=[y1,（y（2:

end）-y（1:

end-1））./h,yend];%f[xn-1,xn]

f2=[f1（2:

end）-f1（1:

end-1）]./[h

（1）,h（1:

end-1）+h（2:

end）,h（end）];%f[xn-1,xn,xn+1]

A=2.*eye（n+1）;

A（2:

end,1:

end-1）=A（2:

end,1:

end-1）+diag（[miu,1]'）;

A（1:

end-1,2:

end）=A（1:

end-1,2:

end）+diag（[1,lamda]'）;

M=A\（6*f2'）;

Sx=[y（1:

end-1）',（（y（2:

end）-y（1:

end-1））./h）'-（（1/3*M（1:

end-1）+1/6*M（2:

end））.*h'）,1/2*M（1:

end-1）,1/6*（M（2:

end）-M（1:

end-1））./h'];

%%

fori=0:

9

xx=i+0.5;

forj=2:

n+1

ifxx

S（i+1）=Sx（j-1,:

）*[1,xx-x（j-1）,（xx-x（j-1））^2,（xx-x（j-1））^3]';

break

end

end

end

x∈[0,1]时；S（x）=2.51+0.8x-0.0014861x2-0.00851395x3

x∈[1,2]时；S（x）=3.3+0.771486（x-1）-0.027028（x-1）2-0.00445799（x-1）3

x∈[2,3]时；S（x）=4.04+0.704056（x-2）-0.0404019（x-2）2-0.0036543（x-2）3

x∈[3,4]时；S（x）=4.7+0.612289（x-3）-0.0513648（x-3）2-0.0409245（x-3）3

x∈[4,5]时；S（x）=5.22+0.386786（x-4）-0.174138（x-4）2+0.107352（x-4）3

x∈[5,6]时；S（x）=5.54+0.360567（x-5）+0.147919（x-5）2-0.268485（x-5）3

x∈[6,7]时；S（x）=5.78-0.149051（x-6）-0.657537（x-6）2+0.426588（x-6）3

x∈[7,8]时；S（x）=5.4-0.184361（x-7）+0.622227（x-7）2-0.267865（x-7）3

x∈[8,9]时；S（x）=5.57+0.256496（x-8）-0.181369（x-8）2+0.0548728（x-8）3

x∈[9,10]时；S（x）=5.7+0.058376（x-9）-0.0167508（x-9）2+0.0583752（x-9）3

S（0.5）=2.90856S（1.5）=3.67843S（2.5）=4.38147

S（3.5）=4.98819S（4.5）=5.38328S（5.5）=5.7237

S（6.5）=5.59441S（7.5）=5.42989S（8.5）=5.65976

S（9.5）=5.7323

习题五重积分的计算

23（上机题）重积分的计算

题目：

给定积分

。

取初始步长h和k，及精度

。

应用复化梯形公式，采用逐次二分步长的方法，编制计算I（f）的通用程序。

计算至相邻两次近似值之差的绝对值不超过

为止。

1）用所编程序计算积分

，取

。

程序：

clc

clear

%%example

f=inline（'tan（x.^2+y.^2）','x','y'）;

a=0;

b=pi/3;

c=0;

d=pi/6;

%%define

error=1;

k=1;

n=1;

while（error>0.5e-5）

[x,y]=meshgrid（linspace（c,d,2^k+1）,linspace（a,b,2^k+1））;

h=（b-a）/2^k;

l=（d-c）/2^k;

z=f（x,y）;

z1=z（1:

end-1,1:

end-1）;

z2=z（1:

end-1,2:

end）;

z3=z（2:

end,1:

end-1）;

z4=z（2:

end,2:

end）;

t（k）=h*l/4*（sum（sum（z1））+sum（sum（z2））+sum（sum（z3））+sum（sum（z4）））;

%%extrapolation

if（k>=2）

T（1,k-1）=4/3*t（k）-1/3*t（k-1）;%T

（1）

error=min（error,abs（t（k）-t（k-1）））;

if（k>=3）

T（2,k-2）=16/15*T（1,k-1）-1/15*T（1,k-2）;%T

（2）

error=min（error,abs（T（1,k-1）-T（1,k-2）））;

if（k>=4）

T（3,k-3）=64/63*T（1,k-2）-1/63*T（1,k-3）;%T（3）

error=min（error,abs（T（2,k-2）-T（2,k-3）））;

if（k>=5）

error=min（error,abs（T（3,k-3）-T（3,k-4）））;

end

end

end

end

k=k+1;

end

计算结果：

T（f）T

（1）（f）T

（2）（f）T（3）（f）

1

0.5197965

0.344032

0.337393

0.337709

2

0.3879734

0.337808

0.336592

0.33665

3

0.3503495

0.336668

0.336524

0.336531

4

0.3400887

0.336533

0.336521

5

0.3374218

0.336521

6

0.3367464

I（f）=0.33652

二分6次

习题6

23．（上机题）常微分方程初值问题数值解

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

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

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

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

（5）对于初值问题

取步长

应用

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

作比较；

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

程序：

clc

clear

%%Originalquestion

f=inline（'-x*x*y*y','x','y'）;

y0=3;

h=0.1;

xstr=0;

xend=1.5;

x=xstr:

h:

xend;

yx=3./（1+x.^3）;

n=size（x,2）;

%%RK4method

RK4y

（1）=y0;

fori=1:

n-1

k1=f（x（i）,RK4y（i））;

k2=f（x（i）+h/2,RK4y（i）+h/2*k1）;

k3=f（x（i）+h/2,RK4y（i）+h/2*k2）;

k4=f（x（i）+h,RK4y（i）+h*k3）;

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

end

%%AB4method

AB4y（1:

4）=RK4y（1:

4）;

fori=4:

n-1

AB4y（i+1）=AB4y（i）+h/24*（55*f（x（i）,AB4y（i））-59*f（x（i-1）,AB4y（i-1））+37*f（x（i-2）,AB4y（i-2））-9*f（x（i-3）,AB4y（i-3）））;

end

%%AB4-AM4predictivemethod

BM4y（1:

4）=RK4y（1:

4）;

fori=4:

n-1

yp（i+1）=BM4y（i）+h/24*（55*f（x（i）,BM4y（i））-59*f（x（i-1）,BM4y（i-1））+37*f（x（i-2）,BM4y（i-2））-9*f（x（i-3）,BM4y（i-3）））;

BM4y（i+1）=BM4y（i）+h/24*（9*f（x（i+1）,yp（i+1））+19*f（x（i）,BM4y（i））-5*f（x（i-1）,BM4y（i-1））+f（x（i-2）,BM4y（i-2）））;

end

%%ImprovedAB4-AM4predictivemethod

imprBM4y（1:

4）=RK4y（1:

4）;

fori=4:

n-1

yP（i+1）=imprBM4y（i）+h/24*（55*f（x（i）,imprBM4y（i））-59*f（x（i-1）,imprBM4y（i-1））+37*f（x（i-2）,imprBM4y（i-2））-9*f（x（i-3）,imprBM4y（i-3）））;

yc（i+1）=imprBM4y（i）+h/24*（9*f（x（i+1）,yP（i+1））+19*f（x（i）,imprBM4y（i））-5*f（x（i-1）,imprBM4y（i-1））+f（x（i-2）,imprBM4y（i-2）））;

imprBM4y（i+1）=251/270*yc（i+1）+19/270*yP（i+1）;

end

%%Error

error（1:

4,1:

n）=abs（[yx-RK4y;yx-AB4y;yx-BM4y;yx-imprBM4y]）;

计算结果：

kx（k）y（x）RK4方法误差AB4方法误差AB4—AM4误差带改进AB4—AM4误差

1

0

3

3

0

3

0

3

0

3

0

2

0.1

2.997003

2.997003

1.87E-07

2.997003

1.87E-07

2.997003

1.87E-07

2.997003

1.87E-07

3

0.2

2.97619

2.97619

3.92E-07

2.97619

3.92E-07

2.97619

3.92E-07

2.97619

3.92E-07

4

0.3

2.92113

2.921129

7.58E-07

2.921129

7.58E-07

2.921129

7.58E-07

2.921129

7.58E-07

5

0.4

2.819549

2.819547

1.61E-06

2.818389

0.00116

2.819678

0.00013

2.819588

3.88E-05

6

0.5

2.666667

2.666663

3.18E-06

2.664672

0.001994

2.666876

0.000209

2.666713

4.62E-05

7

0.6

2.467105

2.4671

5.01E-06

2.465203

0.001903

2.467252

0.000147

2.467097

8.23E-06

8

0.7

2.233805

2.233799

5.77E-06

2.233079

0.000726

2.233731

7.35E-05

2.233682

0.000122

9

0.8

1.984127

1.984123

4.13E-06

1.984951

0.000824

1.983787

0.00034

1.983885

0.000242

10

0.9

1.735107

1.735107

1.16E-07

1.737043

0.001936

1.734607

0.0005

1.734808

0.000299

11

1

1.5

1.500006

5.81E-06

1.502195

0.002195

1.499516

0.000484

1.499732

0.000268

12

1.1

1.287001

1.287013

1.13E-05

1.288763

0.001762

1.286657

0.000344

1.286821

0.000181

13

1.2

1.099707

1.099722

1.54E-05

1.100724

0.001017

1.099533

0.000174

1.099622

8.50E-05

14

1.3

0.93838

0.938397

1.77E-05

0.93871

0.000331

0.938343

3.72E-05

0.938367

1.24E-05

15

1.4

0.801282

0.8013

1.84E-05

0.801135

0.000147

0.801327

4.53E-05

0.801311

2.93E-05

16

1.5

0.685714

0.685732

1.78E-05

0.685335

0.00038

0.685796

8.18E-05

0.68576

4.62E-05

结论：

带改进的AB4-AM4预测校正方法比AB4-AM4预测校正方法精度更高，AB4方法精度最低，RK4方法的精度最高