东南大学数值分析上机题作业MATLAB版.docx
《东南大学数值分析上机题作业MATLAB版.docx》由会员分享，可在线阅读，更多相关《东南大学数值分析上机题作业MATLAB版.docx（27页珍藏版）》请在冰豆网上搜索。
东南大学数值分析上机题作业MATLAB版
11
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（'========================================================'）;
程序
>>P20T17
请输入N值:
10^2
精确值为：
从大到小的顺序累加得SN=
从小到大的顺序累加得SN=
============================================================
>>P20T17
请输入N值:
10^4
精确值为：
从大到小的顺序累加得SN=
运行结果
从小到大的顺序累加得SN=
============================================================
>>P20T17
请输入N值:
10^6
精确值为：
从大到小的顺序累加得SN=
从小到大的顺序累加得SN=
============================================================
结果分析
按从大到小的顺序，有效位数分别为：
6，4，3。
按从小到大的顺序，有效位数分别为：
5，6，6。
可以看出，不同的算法造成的误差限是不同的，好的算法可以让结果更加精确。
当采用从大到小的顺序累加的算法时，误差限随着N的增大而增大，可见在累加的过程中，误差在放大，造成结果的误差较大。
因此，采取从小到大的顺序累加得到的结果更加精确。
2
题目
（1）给定初值
及容许误差
，编制牛顿法解方程f（x）=0的通用程序。
（2）给定方程
易知其有三个根
由牛顿方法的局部收敛性可知存在
当
时，Newton迭代序列收敛于根x2*。
试确定尽可能大的
。
试取若干初始值，观察当
时Newton序列的收敛性以及收敛于哪一个根。
（3）通过本上机题，你明白了什么
程序
函数m文件：
functionFu=fu（x）
Fu=x^3/3-x;
end
函数m文件：
functionFu=dfu（x）
Fu=x^2-1;
end
用Newton法求根的通用程序
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）;
寻找最大δ值的程序：
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）;
运行结果
（1）寻找最大的
值。
算法为：
将初值x0在从0开始不断累加搜索精度eps，带入Newton迭代公式，直到求得的根不再收敛于0为止，此时的x0值即为最大的sigma值。
运行，得到在不同的搜索精度下的最大sigma值。
>>Find
请输入搜索精度：
10^-6
请输入容许误差：
10^-6
最大的sigma值为：
>>Find
请输入搜索精度：
10^-4
请输入容许误差：
10^-6
最大的sigma值为：
>>Find
请输入搜索精度：
10^-2
请输入容许误差：
10^-6
最大的sigma值为：
（2）运行
在
内取初值，运行结果如下：
X0
Xk
-1000
-500
-100
-10
-5
可见，在
区间内取初值，Newton序列收敛，且收敛于根
。
在
内取初值，运行结果如下：
X0
Xk
可见，在
内取初值，Newton序列收敛，且收敛于根
。
在
内内取初值，运行结果如下：
X0
Xk
可见，在
内取初值，Newton序列收敛，且收敛于根0。
在
内取初值，运行结果如下：
X0
Xk
可见，在
内取初值，Newton序列收敛，且收敛于根
在
内取初值，运行结果如下：
X0
Xk
5
10
100
500
1000
可见，在
内取初值，Newton序列收敛，且收敛于根
3
题目
对于某电路的分析，归结为求解线性方程组RI=V，其中
（1）编制解n阶线性方程组
的列主元高斯消去法的通用程序；
（2）用所编程序线性方程组
，并打印出解向量，保留5位有效数字；
（3）本题编程之中，你提高了哪些编程能力
程序
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）;
运行结果
运行程序，输入系数矩阵和方程组右端列向量。
运行过程与结果如下图所示：
>>P126T39
请输入线性方程组阶数:
n=4
请输入系数矩阵：
A=
[00;0;0;00]
请输入线性方程组右端向量：
b=
[]
方程组的解为:
2495
>>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]
方程组的解为:
可看出，算得的该线性方程组的解向量为：
[]
4
题目
（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
端点条件为
，用所编程序求车门的3次样条插值函数S（x）,并打印出S（i+，i=0,1，…，9。
程序
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）+;
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+'）
disp（'ix（i+S（i+'）;
fori=1:
n-1
fprintf（'%d%.5f%.5f\n',i,xn（i）+,S（i））;
end
运行结果
>>P195T37
请输入节点数：
n=11
请输入节点坐标：
[012345678910]
请输入节点处函数值：
[]
请输入左边界条件：
y’（x0）=
请输入右边界条件：
y’（xn）=
S（x）=
*x+*x^2-*x^3+（0======================================================================
*x-*x^2-*x^3+（1======================================================================
*x-*x^2-*x^3+（2======================================================================
*x^2-*x-*x^3+（3======================================================================
*x-*x^2+*x^3-（4======================================================================
*x^2-*x-*x^3+（5======================================================================
*x-*x^2+*x^3-（6======================================================================
*x^2-*x-*x^3+（7======================================================================
*x-*x^2+*x^3-（8======================================================================
*x-*x^2+*x^3-（9======================================================================
S（i+
ix（i+S（i+
1
2
3
4
5
6
7
8
9
10
5
题目
用Romberg求积法计算积分
的近似值，要求误差不超过
。
程序
%被积函数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
运行结果
>>Romberg
请输入积分下限：
a=-1
请输入积分上限：
b=1
请输入允许精度：
eps=*10^-7
该积分的值为：
结果分析
手动化简该定积分并最终求得的值为：
，误差限为：
，可见，程序完成了计算要求。
6
题目
常微分方程初值问题数值解
（1）编制RK4方法的通用程序；
（2）编制AB4方法的通用程序（由RK4提供初值）；
（3）编制AB4-AM4预测校正方法通用程序（由RK4提供初值）；
（4）编制带改进的AB4-AM4预测校正方法通用程序（由RK4提供初值）；
（5）对于初值问题
取步长h=，应用
（1）-（4）中的四种方法进行计算，并将计算结果和精确解
作比较；
（6）通过本上机题，你能得到哪些结论
程序
%f（x，y）函数m文件：
functionFXY=fxy（x,y）
FXY=-x*x*y*y;
end
%精确解y（x）函数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）+*h,y（i）+*h*k1）;
k3=fxy（x（i）+*h,y（i）+*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）+*h,y（i）+*h*k1）;
k3=fxy（x（i）+*h,y（i）+*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*fxy（x（i-2）,y（i-2））-9*fxy（x（i-3）,y（i-3）））;
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-AM4预测校正法通用程序
functionAB4AM4（）
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）+*h,y（i）+*h*k1）;
k3=fxy（x（i）+*h,y（i）+*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
yp（i+1）=y（i）+h/24*（55*fxy（x（i）,y（i））-59*fxy（x（i-1）,y（i-1））+37*fxy（x（i-2）,y（i-2））-9*fxy（x（i-3）,y（i-3）））;
y（i+1）=y（i）+h/24*（9*fxy（x（i+1）,yp（i+1））+19*fxy（x（i）,y（i））-5*fxy（x（i-1）,y（i-1））+fxy（x（i-2）,y（i-2）））;
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））;
d
升级会员 