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

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

2015.1.9

上机作业题报告
USER

1.Chapter 1
1.1题目
设SN=j=2N1j2-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);
for i=2:N
Snl2s=Snl2s+1/(i*i-1);
end
for i=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.Chapter 2
2.1题目
（1）给定初值及容许误差，编制牛顿法解方程f(x)=0的通用程序。
（2）给定方程,易知其有三个根
由牛顿方法的局部收敛性可知存在当时，Newton迭代序列收敛于根x2*。试确定尽可能大的δ。
试取若干初始值，观察当时Newton序列的收敛性以及收敛于哪一个根。
（3）通过本上机题，你明白了什么？
2.2程序
f(x)函数m文件：fu.m
function Fu=fu(x)
Fu=x^3/3-x;
end

f'(x)函数m文件：dfu.m
function Fu=dfu(x)
Fu=x^2-1;
end

用Newton法求根的通用程序Newton.m
clear;
x0=input('请输入初值x0:');
ep=input('请输入容许误差:');
flag=1;
while flag==1
x1=x0-fu(x0)/dfu(x0);
if abs(x1-x0)<ep
flag=0;
end

x0=x1;
end
fprintf('方程的一个近似解为：%f\n',x0);

寻找最大δ值的程序：Find.m
clear
eps=input('请输入搜索精度：');
ep=input('请输入容许误差：');
flag=1;
k=0;
x0=0;
while flag==1
sigma=k*eps;
x0=sigma;
k=k+1;
m=0;
flag1=1;
while flag1==1 && m<=10^3
x1=x0-fu(x0)/dfu(x0);
if abs(x1-x0)<ep
flag1=0;
end
m=m+1;
x0=x1;
end
if flag1==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
在（-∞，-1）内取初值，运行结果如下：
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
可见，在（-∞，-1）区间内取初值，Newton序列收敛，且收敛于根-3。
在（-1，-δ）内取初值，运行结果如下：
X0
Xk
-0.95
1.732051
-0.85
1.732051
-0.8
1.732051
-0.774598
1.732051
可见，在（-1，-δ）内取初值，Newton序列收敛，且收敛于根3。
在（-δ，δ）内内取初值，运行结果如下：
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。
在（δ，1）内取初值，运行结果如下：
X0
Xk
0.774598
-1.732051
0.8
-1.732051
0.85
-1.732051
0.95
-1.732051
可见，在（δ，1）内取初值，Newton序列收敛，且收敛于根-3
在（1，+∞）内取初值，运行结果如下：
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
可见，在（1，+∞）内取初值，Newton序列收敛，且收敛于根3
3.Chapter 3
3.1题目
对于某电路的分析，归结为求解线性方程组RI=V，其中

（1）编制解n阶线性方程组Ax=b的列主元高斯消去法的通用程序；
（2）用所编程序线性方程组RI=V，并打印出解向量，保留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];
for i=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;
for k=i+1:n
if C(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);
for i=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.01 90.86 0 0;90.86 98.81 -67.59 0;0 -67.59 132.01 46.26 ;0 0 46.26 177.17]
请输入线性方程组右端向量：b=
[-33.254 49.79 28.067 -7.324]
方程组的解为:
-2957.4
4426.6
2495
-651.49
>> P126T39
请输入线性方程组阶数: n=9
请输入系数矩阵：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]
方程组的解为:
-0.28923
0.34544
-0.71281
-0.22061
-0.4304
0.15431
-0.057823
0.20105
0.29023
可看出，算得的该线性方程组的解向量为：
[-0.28923 0.34544 -0.71281 -0.22061 -0.4304 0.15431 -0.057823 0.20105 0.29023]
4.Chapter 4
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
端点条件为y0'=0.8， y10'=0.2，用所编程序求车门的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);
for i=1:n-1
h(i)=xn(i+1)-xn(i);
f1(i)=(yn(i+1)-yn(i))/h(i);
end
for i=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);
for i=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;
for i=1:n
A(i,i)=2;
end
for i=2:n-1
A(i,i-1)=miu(i-1);
A(i,i+1)=lamda(i-1);
end
M=A\d;
%====================回代求插值函数====================%
syms x;
for i=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);
for i=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)=');
for i=1:n-1
format short;
fprintf(' %s (%d<x<%d)\n',char(Sx(i)),xn(i),xn(i+1));
disp('======================================================================');
end
disp('S(i+0.5)')
disp(' i x(i+0.5) S(i+0.5)');
for i=1:n-1
fprintf(' %d %.5f %.5f\n',i,xn(i)+0.5,S(i));
end
4.3运行结果
>> P195T37
请输入节点数：n=11
请输入节点坐标：[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]
请输入左边界条件：y’（x0）=0.8
请输入右边界条件：y’（xn）=0.2
S(x)=
0.79*x + 0.0158344*x^2 - 0.0158344*x^3 + 2.51 (0<x<1)
======================================================================
0.830013*x - 0.0241785*x^2 - 0.00249676*x^3 + 2.49666 (1<x<2)
======================================================================
0.809832*x - 0.0140879*x^2 - 0.00417854*x^3 + 2.51012 (2<x<3)
======================================================================
0.315407*x^2 - 0.178653*x - 0.0407891*x^3 + 3.4986 (3<x<4)
======================================================================
6.9313*x - 1.46208*x^2 + 0.107335*x^3 - 5.98133 (4<x<5)
======================================================================
4.1762*x^2 - 21.2601*x - 0.26855*x^3 + 41.0043 (5<x<6)
======================================================================
53.8449*x - 8.3413*x^2 + 0.426866*x^3 - 109.206 (6<x<7)
======================================================================
6.27011*x^2 - 48.435*x - 0.268915*x^3 + 129.447 (7<x<8)
======================================================================
14.4854*x - 1.59494*x^2 + 0.0587951*x^3 - 38.3403 (8<x<9)
======================================================================
13.2458*x - 1.45831*x^2 + 0.053735*x^3 - 34.5615 (9<x<10)
======================================================================
S(i+0.5)
i x(i+0.5) S(i+0.5)
1 0.50000 2.90698
2 1.50000 3.67885
3 2.50000 4.38136
4 3.50000 4.98822
5 4.50000 5.38326
6 5.50000 5.72372
7 6.50000 5.59435
8 7.50000 5.43012
9 8.50000 5.65892
10 9.50000 5.73172
5.Chapter 5
5.1题目
用Romberg求积法计算积分
-1111+100x2dx
的近似值，要求误差不超过0.5×10-7。
5.2程序
%被积函数m文件：fx.m
function Fx=fx(x)
Fx=1/(1+100*x*x);
end

%Romberg求积法计算积分的通用程序
function Romberg()
clear;
a=input('请输入积分下限：a=');
b=input('请输入积分上限：b=');
eps=input('请输入允许精度：eps=');
%========计算Tn========%
function Tn=T(n)
Tn=0;
h=(b-a)/n;
x=zeros(1,n+1);
for k=1:n+1
x(k)=a+(k-1)*h;
end
for j=1:n
Tn=Tn+h*(fx(x(j))+fx(x(j+1)))/2;
end
end
%========计算Sn========%
function Sn=S(n)
Sn=4/3*T(2*n)-1/3*T(n);
end
%========计算Cn========%
function Cn=C(n)
Cn=16/15*S(2*n)-1/15*S(n);
end
%========计算Rn========%
function Rn=R(n)
Rn=64/63*C(2*n)-1/63*C(n);
end
%========计算满足允许精度的Rn，并打印输出========%
i=1;
flag=1;
while flag==1
if abs(R(2^i)-R(2^(i-1)))/255<eps
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，误差限为：3.486×10-8，可见，程序完成了计算要求。
6.Chapter 6
6.1题目
常微分方程初值问题数值解
（1）编制RK4方法的通用程序；
（2）编制AB4方法的通用程序（由RK4提供初值）；
（3）编制AB4-AM4预测校正方法通用程序（由RK4提供初值）；
（4）编制带改进的AB4-AM4预测校正方法通用程序（由RK4提供初值）；
（5）对于初值问题
y'=-x2y2y0=3
取步长h=0.1，应用（1）-（4）中的四种方法进行计算，并将计算结果和精确解yx=3/(1+x3)作比较；
（6）通过本上机题，你能得到哪些结论？
6.2程序
%f（x，y）函数m文件：fxy.m
function FXY=fxy(x,y)
FXY=-x*x*y*y;
end

%精确解y（x）函数m文件：fx.m
function FX=fx(x)
FX=3/(1+x*x*x);
end

%RK4法通用程序
function RK4()
clear;
x(1)=input('请输入初始x值：x0=');
y(1)=input('请输入初值条件：y（x0）=');
N=input('请输入计算步长：N=');
h=input('请输入步长：h=');
for i=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('i xi yi y(xi) y(xi)-yi');
disp('-------------------------------------------------------------');
for i=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法通用程序
function AB4()
clear;
x(1)=input('请输入初始x值：x0=');
y(1)=input('请输入初值条件：y（x0）=');
N=input('请输入计算步长：N=');
h=input('请输入步长：h=');
for i=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
for i=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('i xi yi y(xi) y(xi)-yi');
disp('-------------------------------------------------------------');
for i=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预测校正法通用程序
function AB4AM4()
clear;
x(1)=input('请输入初始x值：x0=');
y(1)=input('请输入初值条件：y（x0）=');
N=input('请输入计算步长：N=');
h=input('请输入步长：h=');
for i=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
for i=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('i xi yi y(xi) y(xi)-yi');
disp('-------------------------------------------------------------');
for i=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预测校正方法
function AB4AM4plus()
clear;
x(1)=input('请输入初始x值：x0=');
y(1)=input('请输入初值条件：y（x0）=');
N=input('请输入计算步长：N=');
h=input('请输入步长：h=');
for i=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
for i=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)));
yc(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)));
y(i+1)=251/270*yc(i+1)+19/270*yp(i+1);
end
disp('i xi yi y(xi) y(xi)-yi');
disp('-------------------------------------------------------------');
for i=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
6.3运行结果
（1）RK4法
>> RK4
请输入初始x值：x0=0
请输入初值条件：y（x0）=3
请输入计算步数：N=15
请输入步长：h=0.1
i xi yi y(xi) y(xi)-yi
-------------------------------------------------------------
1 0.000000 3.000000 3.000000 0.000000
-------------------------------------------------------------
2 0.100000 2.997003 2.997003 0.000000
-------------------------------------------------------------
3 0.200000 2.976190 2.976190 0.000000
-------------------------------------------------------------
4 0.300000 2.921129 2.921130 0.000001
-------------------------------------------------------------
5 0.400000 2.819547 2.819549 0.000002
-------------------------------------------------------------
6 0.500000 2.666663 2.666667 0.000003
-------------------------------------------------------------
7 0.600000 2.467100 2.467105 0.000005
-------------------------------------------------------------
8 0.700000 2.233799 2.233805 0.000006
-------------------------------------------------------------
9 0.800000 1.984123 1.984127 0.000004
-------------------------------------------------------------
10 0.900000 1.735107 1.735107 -0.000000
-------------------------------------------------------------
11