东南大学《数值分析》上机题.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
ifxxS=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
ifxxS（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方法的精度最高
升级会员 