东南大学数值分析上机题作业MATLAB版.docx
《东南大学数值分析上机题作业MATLAB版.docx》由会员分享，可在线阅读，更多相关《东南大学数值分析上机题作业MATLAB版.docx（27页珍藏版）》请在冰豆网上搜索。
东南大学数值分析上机题作业MATLAB版
1.Chapter1
1.1题目
设
，其精确值为
。
（1）编制按从大到小的顺序
，计算SN的通用程序。
（2）编制按从小到大的顺序
，计算SN的通用程序。
（3）按两种顺序分别计算
并指出有效位数。
（编制程序时用单精度）
（4）通过本次上机题，你明白了什么？
clear;
N=input（'请输入N值:
'）;
Ac=single（（3/2-1/N-1/（N+1））/2）;
Snl2s=single（0）;
Sns2l=single（0）;
fori=2:
N
Snl2s=Snl2s+1/（i*i-1）;
end
fori=N:
-1:
2
Sns2l=Sns2l+1/（i*i-1）;
end
fprintf（'精确值为：
%f\n',Ac）;
fprintf（'从大到小的顺序累加得SN=%f\n',Snl2s）;
fprintf（'从小到大的顺序累加得SN=%f\n',Sns2l）;
disp（'========================================================'）;
1.2程序
>>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（%ddisp（'======================================================================'）;
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）））/255flag=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.4860747，误差限为：
，可见，程序完成了计算要求。
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（'--------------------------------
升级会员 