数值实验大作业.docx
- 文档编号:24541134
- 上传时间:2023-05-28
- 格式:DOCX
- 页数:29
- 大小:715.47KB
数值实验大作业.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
whilep a=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 whilep a=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');%精确解图像
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 数值 实验 作业