北理工数值分析大作业.docx
《北理工数值分析大作业.docx》由会员分享,可在线阅读,更多相关《北理工数值分析大作业.docx(24页珍藏版)》请在冰豆网上搜索。
北理工数值分析大作业
数值分析上机作业
第1章
1.1计算积分
,n=9。
(要求计算结果具有6位有效数字)
程序:
n=1:
19;
I=zeros(1,19);
I(19)=1/2*((exp(-1)/20)+(1/20));
I(18)=1/2*((exp(-1)/19)+(1/19));
fori=2:
10
I(19-i)=1/(20-i)*(1-I(20-i));
end
formatlong
disp(I(1:
19))
结果截图及分析:
在MATLAB中运行以上代码,得到结果如下图所示:
当计算到数列的第10项时,所得的结果即为
n=9时的准确积分值。
取6位有效数字可得
.
1.2分别将区间[-10.10]分为100,200,400等份,利用mesh或surf命令画出二元函数
z=
的三维图形。
程序:
>>x=-10:
0.1:
10;
y=-10:
0.1:
10;
[X,Y]=meshgrid(x,y);
Z=exp(-abs(X))+cos(X+Y)+1./(X.^2+Y.^2+1);
subplot(2,2,1);
mesh(X,Y,Z);
title('步长0.1')
>>x=-10:
0.2:
10;
y=-10:
0.2:
10;
[X,Y]=meshgrid(x,y);
Z=exp(-abs(X))+cos(X+Y)+1./(X.^2+Y.^2+1);
subplot(2,2,1);
mesh(X,Y,Z);
title('步长0.2')
>>x=-10:
0.05:
10;
y=-10:
0.05:
10;
[X,Y]=meshgrid(x,y);
Z=exp(-abs(X))+cos(X+Y)+1./(X.^2+Y.^2+1);
subplot(2,2,1);
mesh(X,Y,Z);
title('步长0.05')
结果截图及分析:
由图可知,步长越小时,绘得的图形越精确。
第2章
试用MATLAB编程实现追赶法求三对角方程组的算法,并考虑梯形电路电阻问题:
电路中的电流
满足下列线性方程组:
设
,求各段电路的电流量。
处理思路:
观察该方程的系数矩阵可知,它是一个三对角矩阵,故可运用追赶法对其进行求解。
程序:
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
结果截图及分析:
在MATLAB中运行以上代码,得到结果如下图所示:
图中8个值依次为
的数值。
第3章
试分别用
(1)Jacobi迭代法;
(2)Gauss-Seidel解线性方程组
迭代初始向量取
.
3.1Jacobi迭代法
程序:
>>A=[101234;
19-12-3;
2-173-5;
32312-1;
4-3-5-115];
b=[12;-27;14;-17;12];
x0=[0;0;0;0;0];
D=diag(diag(A));
I=eye(5);
L=-tril(A,-1);
B=I-D\A;
g=D\b;
y=B*x0+g;
n=1;
whilenorm(y-x0)>=1.0e-6
x0=y;
y=B*x0+g;
n=n+1;
end
fprintf('%8.6f\n',y);
n
结果截图及分析:
得到此结果时迭代次数为67次,达到精度要求。
3.2Gauss-Seidel迭代法:
程序:
>>A=[101234;
19-12-3;
2-173-5;
32312-1;
4-3-5-115];
b=[12;-27;14;-17;12];
x0=[0;0;0;0;0];
D=diag(diag(A));
U=-triu(A,1);
L=-tril(A,-1);
M=(D-L)\U;
g=(D-L)\b;
y=M*x0+g;
n=1;
whilenorm(y-x0)>=1.0e-6
x0=y;
y=M*x0+g;
n=n+1;
end
fprintf('%8.6f\n',y);
结果截图及分析:
Gauss-Seidel迭代法只需要迭代38次即可满足精度要求。
第4章
设A=
,取
先用幂法迭代3次,得到A的按模最大特征值的近似值,取
为其整数部分,再用反幂法计算A的按模最大特征值的更精确的近似值,要求误差小于
.
程序:
A=[126-6;
6162;
-6216];
x0=[1;1;1];
y=x0;b=max(abs(x0));k=1;
while(k<4)
x=A*y;b=max(abs(x));y=x./b;
k=k+1;
fprintf('eig1equals%6.4f\n',b);
end
>>bb0=fix(b);
I=eye(3,3);
x0=[1;1;1];
y=x0;l=0;bb=max(abs(x0));k=1;
while(abs(bb-l)>=1.0e-10)
l=bb;
x=(A-bb0*I)\y;bb=max(abs(x));y=x./bb;
eig=l+b;
>>fprintf('eig2(%d)equals%12.10f\n',k,eig);
k=k+1;
end
实验截图及分析:
由图可知,由幂法3次迭代后得到的特征值为19.4,而由反幂法得到的特征值为20.3999999999.误差小于
第5章
试编写MATLAB函数实现Newton插值,要求能输出插值多项式。
对函数f(x)=
在区间[-5,5]上实现10次多项式插值。
要求:
(1)输出插值多项式。
(2)在区间[-5,5]内均匀插入99个节点,计算这些节点上函数f(x)的近似值,并在同一张图上画出原函数和插值多项式的图形。
(3)观察龙格现象,计算插值函数在各节点处的误差,并画出误差图。
5.1输出插值多项式
程序:
x=-5:
1:
5;
y=1./(1+4*(x.^2));
newpoly(x,y)
function[c,d]=newpoly(x,y)
n=length(x);
d=zeros(n,n);
d(:
1)=y';
forj=2:
n
fork=j:
n
d(k,j)=(d(k,j-1)-d(k-1,j-1))/(x(k)-x(k-j+1));
end
end
c=d(n,n);
fork=(n-1):
-1:
1
c=conv(c,poly(x(k)));
m=length(c);
c(m)=c(m)+d(k,k);
end
end
结果及分析:
ans=
Columns1through2
-0.000049595763049
Columns3through4
0.0027401659084830.000000000000000
Columns5through6
.0514********
Columns7through8
0.3920149852823120.000000000000000
Columns9through10
-1.1432840483510250.000000000000001
Column11
1.000000000000000
10次插值多项式由高到低系数为Columns1至Column11
5.2原函数与插值多项式的图形
程序:
x=-5:
1:
5;
y=1./(1+4*(x.^2));
n=newpoly(x,y);
x0=-5:
0.1:
5;
y0=1./(1+4*(x0.^2));
vn=polyval(n,x0);
plot(x0,vn,'-r',x0,y0,'--b');
xlabel('x');ylabel('y');
实验结果截图:
原函数与插值多项式的图形如上图所示,蓝色为原函数的图形,红色为插值多项式的图形。
5.3各节点的误差及误差图
程序:
formatlong;
x=-5:
1:
5;
y=1./(1+4*(x.^2));
n=newpoly(x,y);
x0=-5:
0.1:
5;
y0=1./(1+4*(x0.^2));
vn=polyval(n,x0);
plot(x0,y0-vn,'-r');
xlabel('x');ylabel('y');
实验结果截图:
误差图如上图所示。
第6章
炼钢厂出钢时所用的盛钢水的钢包,在使用过程中由于钢液及炉渣对包衬耐火材料的腐蚀,使其容积不断加大。
经试验,钢包的容积与相应的使用次数的数据列表如下:
使用次数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,双曲线转化为一次线性方程,使用最小二乘法求出该一次方程的系数。
程序:
x=[2356791011121416171920];
y=[106.42108.26109.58109.5109.86110109.93110.59110.60110.72110.9110.76111.1111.3];
k1=0;k2=0;k3=0;k4=0;
fori=1:
14
k1=k1+1/x(i);
end
fori=1:
14
k2=k2+1/y(i);
end
fori=1:
14
k3=k3+1/(x(i))^2;
end
fori=1:
14
k4=k4+1/(x(i)*y(i));
end
b=(k1*k2-14*k4)/(k1^2-14*k3)
a=k2/14-k1*b/14
plot(x,y,'r*')
holdon
x=2:
0.01:
20;
y=1./(a+b./x);
plot(x,y)
xlabel('x')
ylabel('y')
gridon
实验结果截图与分析:
即最小二乘法求出拟合函数为:
=0.008973+0.000842
拟合曲线图为:
第7章
考纽螺线的形状象钟表的发条,也称回旋曲线,它在直角坐标系中的参数方程为
曲线关于原点对称,取a=1,参数s的变化范围[-5,5],容许误差限分别是
和
。
选取适当的节点个数,利用数值积分方法计算曲线上点的坐标,并画出曲线的图形。
误差限为
时:
程序:
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-6);
i=i+1;
end
plot(x,y,'r-');
title('误差限-1e-6');
xlabel('x(s)');
ylabel('y(s)');
实验截图:
误差限为
时得到曲线如下图所示:
误差限为
时:
程序:
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-10);
y(i,1)=quad(f2,0,s,1e-10);
i=i+1;
end
plot(x,y,'r-');
title('误差限-1e-10');
xlabel('x(s)');
ylabel('y(s)');
实验结果截图:
第8章
求方程
的非零根。
程序:
x=sym('x');
f=sym('log((513+0.6651*x)/(513-0.6651*x))-x/(1400*0.0981)');
df=diff(f,x);
FX=x-f/df;
Fx=inline(FX);
formatlong;
i=1;
x0=760;
fori=1:
10
disp(x0);
x0=feval(Fx,x0);
end
x0
结果截图及分析:
取初值为760,迭代9次,两个非零根为765.48及-765.48.
第9章
设有常微分方程的初值问题
其精确解为
。
选取步长h使四阶Adams预测-校正算法和经典RK法均稳定,分别用这两种方法求解微分方程,将数值解与精确解进行比较,输出结果。
其中多步法需要的初值由经典RK法提供。
程序:
f=inline('-y+2*cos(x)','x','y');
n=20;
Y=zeros(1,n+1);
Y
(1)=1;
h=0.05*pi;
X=0:
h:
pi;
fori=1:
20
K1=f(X(i),Y(i));
K2=f(X(i)+h/2,Y(i)+(h*K1)/2);
K3=f(X(i)+h/2,Y(i)+(h*K2)/2);
K4=f(X(i)+h,Y(i)+h*K3);
Y(i+1)=Y(i)+h*(K1+2*K2+2*K3+K4)/6;
end;
Y1=cos(X)+sin(X);
disp('步长经典RK法精确解');
disp([X'Y'Y1']);
plot(X,Y1,'k*',X,Y,'-r');
gridon;
title('经典RK法解非线性方程');
legend('精确解','经典RK法');
实验结果截图:
四阶Adams预测-校正算法
程序:
f=inline('-y+2*cos(x)','x','y');
n=20;
Y=zeros(1,n+1);
Y
(1)=1;
h=pi/n;
X=0:
h:
pi;
fori=1:
3
K1=f(X(i),Y(i));
K2=f(X(i)+h/2,Y(i)+(h*K1)/2);
K3=f(X(i)+h/2,Y(i)+(h*K2)/2);
K4=f(X(i)+h,Y(i)+h*K3);
Y(i+1)=Y(i)+h*(K1+2*K2+2*K3+K4)/6;
end;
Y1=Y;
fori=4:
n
f1=f(X(i-3),Y(i-3));
f2=f(X(i-2),Y(i-2));
f3=f(X(i-1),Y(i-1));
f4=f(X(i),Y(i));
Y(i+1)=Y(i)+h*(55*f4-59*f3+37*f2-9*f1)/24;
Yn=f(X(i+1),Y(i+1));
Y1(i+1)=Y(i)+h*(9*Yn+19*f4-5*f3+f2)/24;
end
Y2=cos(X)+sin(X);
disp('步长四阶Adams预测-校正算法精确解');
disp([X'Y1'Y2']);
plot(X,Y2,'k*',X,Y1,'-r');
gridon;
title('四阶Adams预测-校正算法解微分方程');
legend('精确解','四阶Adams预测-校正算法');
实验结果截图: