数值实验大作业.docx
《数值实验大作业.docx》由会员分享，可在线阅读，更多相关《数值实验大作业.docx（29页珍藏版）》请在冰豆网上搜索。
数值实验大作业
第一章
11.设
，计算
。
（要求计算结果保留6位有效数字）
解：
Matlab程序：
clearall;
formatlong;
symsx;%定义符号变量
y=int（x^9*exp（x-1）,x,0,1）;%计算定积分
vpa（y）%将结果转换成小数形式
运行结果如下：
ans=
.0916********
分析：
要使结果保留6位有效数字，那么：
4.分别将区间
分为100,200,400等份，利用mesh或surf命令画出二元函数
的三维图形。
解：
Matlab程序：
clearall;
%100等份
x=[-10:
20/100:
10];
y=x;
[X,Y]=meshgrid（x,y）;%生产网格坐标
Z=exp（-abs（X））+cos（X+Y）+1./（X.^2+Y.^2+1）;
figure
（1）;
subplot（1,2,1）;
mesh（X,Y,Z）;
title（'100等份mesh图'）;
subplot（1,2,2）;
surf（X,Y,Z）;
title（'100等份surf图'）;
%200等份
x=[-10:
20/200:
10];
y=x;
[X,Y]=meshgrid（x,y）;%生产网格坐标
Z=exp（-abs（X））+cos（X+Y）+1./（X.^2+Y.^2+1）;
figure
（2）;
subplot（1,2,1）;
mesh（X,Y,Z）;
title（'200等份mesh图'）;
subplot（1,2,2）;
surf（X,Y,Z）;
title（'200等份surf图'）;
%400等份
x=[-10:
20/400:
10];
y=x;
[X,Y]=meshgrid（x,y）;%生产网格坐标
Z=exp（-abs（X））+cos（X+Y）+1./（X.^2+Y.^2+1）;
figure（3）;
subplot（1,2,1）;
mesh（X,Y,Z）;
title（'400等份mesh图'）;
subplot（1,2,2）;
surf（X,Y,Z）;
title（'400等份surf图'）;
运行结果如下：
分析：
由运行结果可知：
1）网格划分越密，得到的三维图形越丰满越精细，越接近函数真实的三维图形。
2）在Matlab中，mesh和surf命令的相同点是：
都可以绘出某一区间内的完整曲面，它们的调用方法类似；不同点是mesh命令绘制的图形是一个一排排的彩色曲线组成的网格图，而surf命令绘制得到的是着色的曲面图，而不是网格近似图。
第二章
1.试用Matlab软件编程实现矩阵的列主元素三角分解，并求Pascal矩阵
的逆阵
，用左除命令A/E检验你的结果。
解：
Matlab程序：
clearall;
formatshort;
A=[1,1,1,1,1;1,2,3,4,5;1,3,6,10,15;1,4,10,20,35;1,5,15,35,70];
[L,U,P]=lu（A）%列主元素三角分解
B=inv（A）%求A的逆
E=eye（5）;
C=A\E%左除
运行结果如下：
分析：
由运行结果可知，B和C相同，即用Matlab命令inv（）求矩阵逆与用左除单位阵得到的结果相同。
2.试用Matlab软件编程实现追赶法求解三对角方程组的算法，并考虑如下梯形电阻电路问题，电路中的各个电流
满足下列线性方程组：
设
，求各段电路的电流量。
解：
Matlab程序：
clearall;
fori=1:
8
a（i）=-2;b（i）=5;c（i）=-2;d（i）=0;
end
a
（1）=0;b
（1）=2;c（8）=0;d
（1）=220/27;
fori=2:
8
a（i）=a（i）/b（i-1）;
b（i）=b（i）-c（i-1）*a（i）;
d（i）=d（i）-a（i）*d（i-1）;
end
d（8）=d（8）/b（8）;
fori=7:
-1:
1
d（i）=（d（i）-c（i）*d（i+1））/b（i）;
end
fori=1:
8
x（i）=d（i）;
end
x
运行结果如下：
x=
8.1477751669091054.0737010928350302.0364775651784711.0174928201111480.5072544850994000.2506433926373500.1193539964939760.047741598597591
分析：
运行结果矩阵x中8个元素分别对应电流
的数值。
用追赶法求解三对角矩阵，能够大大减少计算量。
第三章
1.试分别用
（1）Jacobi迭代法；
（2）Gauss-Seidel迭代法解线性方程组
迭代初始向量取
.
解：
（1）Jacobi迭代法：
首先，判断本方程组的Jacobi迭代法是否收敛，即判断迭代矩阵的谱半径是否小于1。
Matlab程序：
1）判断是否收敛的Matlab程序：
clearall;
A=[10,1,2,3,4;1,9,-1,2,-3;2,-1,7,3,-5;3,2,3,12,-1;4,-3,-5,-1,15];
D=diag（diag（A））;%对角阵
I=eye（5）;%5阶单位阵
M=I-inv（D）*A;%求迭代矩阵
R=vrho（M）%求迭代矩阵的谱半径
运行结果如下：
R=
0.826786213856694
分析：
迭代矩阵M的谱半径R<1，收敛。
2）Jacobi迭代法Matlab程序：
clearall;
%判断迭代法是否收敛
A=[10,1,2,3,4;1,9,-1,2,-3;2,-1,7,3,-5;3,2,3,12,-1;4,-3,-5,-1,15];
D=diag（diag（A））;%对角阵
I=eye（5）;%5阶单位阵
M=I-inv（D）*A;%求迭代矩阵
R=vrho（M）;%求迭代矩阵的谱半径,判断是否小于1
%应用Jacobi迭代法求解
b=[12;-27;14;-17;12];
x0=[0;0;0;0;0];%初始化
g=inv（D）*b;
x=M*x0+g;%Jacobi迭代公式
k=1;%第一次迭代
whilenorm（x-x0）>=1.0e-6%自定义的误差容限
x0=x;
x=M*x+g;
k=k+1;%迭代次数
end
k
x
运行结果如下：
k=
67
x=
1.000001624057766
-2.000001488229411
2.999997271268256
-1.999999562383379
0.999998051336040
分析：
由运行结果可知，使用Jacobi迭代法进行了67次迭代满足自定义的误差容限的要求，得到矩阵x结果如上。
（2）Gauss-Seidel迭代法
首先，判断本方程组的Gauss-Seidel迭代法是否收敛，即判断迭代矩阵的谱半径是否小于1。
Matlab程序：
1）判断是否收敛的Matlab程序：
clearall;
A=[10,1,2,3,4;1,9,-1,2,-3;2,-1,7,3,-5;3,2,3,12,-1;4,-3,-5,-1,15];
D=diag（diag（A））;%对角阵
L=-tril（A,-1）;%A的严格下三角阵的负
U=-triu（A,1）;%A的严格上三角阵的负
M=inv（D-L）*U;%求迭代矩阵
R=vrho（M）%求迭代矩阵的谱半径
运行结果如下：
R=
0.684558446281366
分析：
迭代矩阵M的谱半径R<1，收敛。
2）Gauss-Seidel迭代法Matlab程序：
clearall;
A=[10,1,2,3,4;1,9,-1,2,-3;2,-1,7,3,-5;3,2,3,12,-1;4,-3,-5,-1,15];
D=diag（diag（A））;%对角阵
L=-tril（A,-1）;%A的严格下三角阵的负
U=-triu（A,1）;%A的严格上三角阵的负
M=inv（D-L）*U;%求迭代矩阵
R=vrho（M）;%判断迭代矩阵的谱半径是否小于1
b=[12;-27;14;-17;12];
x0=[0;0;0;0;0];%初始化
g=inv（D-L）*b;
x=M*x0+g;%迭代公式
k=1;%第一次迭代
whilenorm（x-x0）>=1.0e-6%自定义的误差容限
x0=x;
x=M*x+g;
k=k+1;%迭代次数
end
k
x
运行结果如下：
k=
38
x=
1.000000908757070
-2.000000746037326
2.999998708943986
-1.999999879156684
0.999999186161532
分析：
由运行结果可知，使用Gauss-Seidel迭代法进行了38次迭代满足自定义的误差容限的要求，得到矩阵x结果如上。
（3）方程组的精确解
clearall;
formatshort;
A=[10,1,2,3,4;1,9,-1,2,-3;2,-1,7,3,-5;3,2,3,12,-1;4,-3,-5,-1,15];
b=[12;-27;14;-17;12];
x=inv（A）*b
运行结果如下：
x=
1.0000
-2.0000
3.0000
-2.0000
1.0000
总结分析：
对于本题来说，Gauss-Seidel迭代法和Jacobi迭代法都收敛。
Jacobi迭代法进行了67次迭代满足自定义的误差容限的要求，而使用Gauss-Seidel迭代法进行了38次迭代即满足自定义的误差容限的要求，即是说，对于本题，Gauss-Seidel迭代法比Jacobi迭代法收敛得快。
第四章
2.设
，取
，先用幂法迭代3次，得到A的按模最大特征值的近似值，取
为其整数部分，再用反幂法计算A的按模最大特征值的更精确的近似值，要求误差小于
。
解：
幂法Matlab程序：
clearall;
A=[126-6;6162;-6216];
x0=[1;1;1];%初始化
y=x0;
x=A*y;
k=0;
whilek<3
a=max（x）;
y=x/a;%归一化或标准化
x=A*y;%迭代
k=k+1;%记录迭代次数
fprintf（'k=%d:
themaineigis:
%f\n',k,a）;
end
运行结果如下：
k=1:
themaineigis:
24.000000
k=2:
themaineigis:
20.000000
k=3:
themaineigis:
19.400000
分析：
由运行结果可知，使用幂法经过三次迭代，得到A的按膜最大的特征值的近似值为19.4，取其整数部分，
=19。
反幂法Matlab程序：
clearall;
formatlong
A=[12,6,-6;6,16,2;-6,2,16];
d=20;%开始取的19得到的结果不合理，在20附近寻找更精确的主特征值，加速迭代过程
I=eye（3）;
B=A-d*I;%原点移位
x=[1;1;1];%初始化
y=x;
a=1/d;
b=max（x）;
ep=abs（1/a-1/b）;
whileep>1e-10%误差容限
y=x/a;%归一化或标准化
x=inv（B）*y;%迭代公式
a=max（x）;
ep=abs（1/a-1/b）;%相邻两次迭代的误差
b=a;
end
d=d+1/b
运行结果如下：
d=
21.544003745319511
分析：
由运行结果可知，利用反幂法在20附近迭代搜寻到主特征值为21.544003745319511，满足误差容限要求。
本题的精确解：
求特征值的Matlab程序：
clearall;
A=[12,6,-6;6,16,2;-6,2,16];
[V,D]=eig（A）
运行结果如下：
V=
-0.7473423402953060.000000000000001-0.664439181868389
0.469829451185180-0.707106781186547-0.528450836690636
-0.469829451185180-0.7071067811865480.528450836690635
D=
4.45599625468246800
018.0000000000000040
0021.544003745317532
分析：
利用Matlab库函数eig（）可求得矩阵A的特征值及对应的特征向量，由运行结果可知，矩阵A的主特征值为21.544003745317532，而由原点移位反幂法求得的主特征值为21.544003745319511，具有13位有效数字，达到了很高的精度。
第五章
1．试编写Matlab函数实现Newton插值，要求能输出插值多项式。
对函数
在区间[-5,5]上实现10次多项式插值。
（1）输出插值多项式。
（2）在区间[-5,5]内均匀插入99个节点，计算这些节点上函数
的近似值，并在同一张图上画出原函数和插值多项式的图形。
（3）观察龙格现象，计算插值函数在各节点处的误差，并画出误差图。
解：
（1）
Matlab程序：
clearall;
formatlong
x1=[-5:
5];
y1=1./（1+4.*x1.*x1）;
f=y1;
n=11;
i=1;
whilei<11
whilen>i
f（n）=（f（n-1）-f（n））/（x1（n-i）-x1（n））;
n=n-1;
end
i=i+1;
n=11;
end
a=1;
q=1;
p=1;
symsx;
newtonpoly=0;
whileq<12
whilepa=a*（x-x1（p））;
p=p+1;
end
newtonpoly=newtonpoly+f（q）*a;
a=1;
q=q+1;
p=1;
end
newtonpoly
运行结果如下：
newtonpoly=
（36*x）/6565+（3550298616520539*（x+4）*（x+5））/1152921504606846976+（2689247898264063*（x+3）*（x+4）*（x+5））/1152921504606846976+（180********08661*（x+2）*（x+3）*（x+4）*（x+5））/576460752303423488+（462354082176629*（x+1）*（x+2）*（x+3）*（x+4）*（x+5））/144115188075855872-（5850230976024283*x*（x+1）*（x+2）*（x+3）*（x+4）*（x+5））/1152921504606846976+（6518522480310501*x*（x-1）*（x+1）*（x+2）*（x+3）*（x+4）*（x+5））/230584*********3952-（2258610859405831*x*（x-1）*（x+1）*（x-2）*（x+2）*（x+3）*（x+4）*（x+5））/230584*********3952+（1143600435142193*x*（x-1）*（x+1）*（x-2）*（x+2）*（x-3）*（x+3）*（x+4）*（x+5））/4611686018427387904-（7319042784910035*x*（x-1）*（x+1）*（x-2）*（x+2）*（x-3）*（x+3）*（x-4）*（x+4）*（x+5））/147573952589676412928+49/1313
分析：
由运行结果可知，函数
在区间[-5,5]上的10次多项式插值如上。
（2）
Matlab程序：
clearall
formatlong
x1=[-5:
5];
y1=1./（1+4.*x1.*x1）;
x2=[-5:
0.1:
5];
y2=1./（1+4.*x2.*x2）;
plot（x2,y2,'r'）;%原函数用红色线表示
holdon;%保留原图像，不被覆盖
gridon;%打开网格
f=y1;
k=11;
i=1;
whilei<11
whilek>i
f（k）=（f（k-1）-f（k））/（x1（k-i）-x1（k））;
k=k-1;
end
i=i+1;
k=11;
end
i=1;
whilei<=101
a=1;
q=1;
p=1;
x（i）=-5+0.1*（i-1）;
globalnewtonpoly
newtonpoly=0;
whileq<12
whilepa=a*（x（i）-x1（p））;
p=p+1;
end
newtonpoly=newtonpoly+f（q）*a;
a=1;
q=q+1;
p=1;
end
y（i）=newtonpoly;
i=i+1;
end
plot（x,y,'b'）;%牛顿插值多项式函数用蓝色线表示
运行结果如下：
分析：
由运行结果可知，该图像在区间两端存在很大误差，既有龙格现象。
（3）
Matlab程序：
在前面程序的基础上，再输入：
plot（x,abs（y-y2）,'g'）%误差用绿色线表示
运行结果如下：
分析：
由运行结果可知，该图像在区间中部误差较小，但在两端存在很大误差，既有龙格现象。
第六章
2.炼钢厂出钢时所用的盛钢水的钢包，在使用过程中由于钢液及炉渣对包衬耐火材料的侵蚀，使其容积不断增大，经试验，钢包的容积与相应的使用次数的数据列表如下：
使用次数x
容积y
使用次数x
容积y
2
106.42
11
110.59
3
108.26
12
110.60
5
109.58
14
110.72
6
109.50
16
110.90
7
109.86
17
110.76
9
110.00
19
111.10
10
109.93
20
111.30
选用双曲线
对数据进行拟合，使用最小二乘法求出拟合函数，作出拟合曲线图。
解：
首先，将其线性化：
用Y替代1/y，用X替代1/x，原曲线化为Y=a+bX，双曲线转化为一次线性方程，使用最小二乘法求出该一次方程的系数。
Matlab程序：
clearall
x=[2,3,5,6,7,9,10,11,12,14,16,17,19,20];
y=[106.42,108.26,109.58,109.50,109.86,110.00,109.93,110.59,110.60,110.72,110.90,110.76,111.10,111.30];
i=1;
whilei<15
plot（x（i）,y（i）,'r*'）;%描点
holdon;
i=i+1;
end
x=1./x;
y=1./y;%线性化处理
p=polyfit（x,y,1）;%一次多项式拟合
y=polyval（p,x）;%一次拟合函数在x处的取值
x=1./x;
y=1./y;%转换回来
gridon%打开网格
plot（x,y,'g'）%拟合曲线用绿色表示
运行结果如下：
分析：
利用数值计算方法求解：
线性化处理：
取
则有
列写数据表：
0.5
0.009397
0.090909
0.009042
0.333333
0.009237
0.083333
0.009042
0.2
0.009126
0.071429
0.009032
0.166667
0.009132
0.0625
0.009017
0.142857
0.009102
0.058824
0.009029
0.111111
0.009091
0.052632
0.009001
0.1
0.009097
0.05
0.008985
由正则方程得到：
解得：
而由Matlab得到的拟合系数为：
两种方法得到的结果契合。
第七章
1.考纽螺线的形状象钟表的发条，也称回旋曲线，它在直角坐标系中的参数方程为
曲线关于原点对称，取a=1，参数s的变化范围[-5,5]，容许误差限分别是
和
。
选取适当的节点个数，利用数值积分方法计算曲线上点的坐标，并画出曲线的图形。
解：
Matlab程序：
clearall
x=zeros（101,1）;
y=zeros（101,1）;
f1=inline（'cos（1/2*（t.^2））'）;
f2=inline（'sin（1/2*（t.^2））'）;%定义函数，方便后续使用
i=1;
fors=-5:
0.1:
5
x（i,1）=quad（f1,0,s,1e-6）;
y（i,1）=quad（f2,0,s,1e-10）;%积分上限是变量
i=i+1;
end
plot（x,y,'b'）;
xlabel（'x（s）'）;
ylabel（'y（s）'）;%横纵坐标
运行结果如下：
第九章
1.设有常微分方程的初值问题
其精确解为
。
选取步长h使四阶Adams预测-校正算法和经典RK法均稳定，分别用这两种方法求解微分方程，将数值解与精确解进行比较，输出结果。
其中多步法需要的初值由经典RK法提供。
解：
Matlab程序：
1）经典4阶RK方法：
clearall
%RK4方法
a=0;
b=pi;
n=100;%100等分
h=（b-a）/n;%步长
x=a:
h:
b;
y=zeros（1,n）;
y
（1）=1;%初始条件
%Rk4方法近似求解
%计算
fori=1:
n-1
k1=-y（i）+2*cos（x（i））;
k2=-（y（i）+h/2*k1）+2*cos（x（i）+h/2）;
k3=-（y（i）+h/2*k2）+2*cos（x（i）+h/2）;
k4=-（y（i）+h*k3）+2*cos（x（i）+h）;
y（i+1）=y（i）+h/6*（k1+2*k2+2*k3+k4）;
end
%出图
fori=1:
n
plot（x（i）,y（i）,'ro'）;%标出各节点
holdon;
end
%精确解
y1=cos（x）+sin（x）;
plot（x,y1,'b'）;%精确解图像
升级会员 