北理工数值分析大作业概要Word文档下载推荐.docx
《北理工数值分析大作业概要Word文档下载推荐.docx》由会员分享,可在线阅读,更多相关《北理工数值分析大作业概要Word文档下载推荐.docx(26页珍藏版)》请在冰点文库上搜索。
mesh(X,Y,Z);
title('
步长0.1'
)
0.2:
步长0.2'
x=-10:
0.05:
步长0.05'
由图可知,步长越小时,绘得的图形越精确。
第2章
试用MATLAB编程实现追赶法求三对角方程组的算法,并考虑梯形电路电阻问题:
电路中的电流
满足下列线性方程组:
设
,求各段电路的电流量。
处理思路:
观察该方程的系数矩阵可知,它是一个三对角矩阵,故可运用追赶法对其进行求解。
程序:
fori=1:
8
a(i)=-2;
b(i)=5;
c(i)=-2;
d(i)=0;
a
(1)=0;
b
(1)=2;
c(8)=0;
d
(1)=220/27;
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);
d(8)=d(8)/b(8);
fori=7:
-1:
1
d(i)=(d(i)-c(i)*d(i+1))/b(i);
x(i)=d(i);
x
图中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];
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;
fprintf('
%8.6f\n'
y);
n
得到此结果时迭代次数为67次,达到精度要求。
3.2Gauss-Seidel迭代法:
D=diag(diag(A));
U=-triu(A,1);
M=(D-L)\U;
g=(D-L)\b;
y=M*x0+g;
y=M*x0+g;
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;
eig1equals%6.4f\n'
b);
bb0=fix(b);
I=eye(3,3);
l=0;
bb=max(abs(x0));
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);
实验截图及分析:
由图可知,由幂法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:
fork=j:
d(k,j)=(d(k,j-1)-d(k-1,j-1))/(x(k)-x(k-j+1));
end
c=d(n,n);
fork=(n-1):
c=conv(c,poly(x(k)));
m=length(c);
c(m)=c(m)+d(k,k);
结果及分析:
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原函数与插值多项式的图形
n=newpoly(x,y);
x0=-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;
plot(x0,y0-vn,'
误差图如上图所示。
第6章
炼钢厂出钢时所用的盛钢水的钢包,在使用过程中由于钢液及炉渣对包衬耐火材料的腐蚀,使其容积不断加大。
经试验,钢包的容积与相应的使用次数的数据列表如下:
使用次数x
容积y
使用次数x
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
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;
k1=k1+1/x(i);
k2=k2+1/y(i);
k3=k3+1/(x(i))^2;
k4=k4+1/(x(i)*y(i));
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)
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:
x(i,1)=quad(f1,0,s,1e-6);
y(i,1)=quad(f2,0,s,1e-6);
i=i+1;
r-'
误差限-1e-6'
x(s)'
y(s)'
实验截图:
时得到曲线如下图所示:
x(i,1)=quad(f1,0,s,1e-10);
y(i,1)=quad(f2,0,s,1e-10);
误差限-1e-10'
第8章
求方程
的非零根。
x=sym('
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);
x0=760;
disp(x0);
x0=feval(Fx,x0);
x0
取初值为760,迭代9次,两个非零根为765.48及-765.48.
第9章
设有常微分方程的初值问题
其精确解为
选取步长h使四阶Adams预测-校正算法和经典RK法均稳定,分别用这两种方法求解微分方程,将数值解与精确解进行比较,输出结果。
其中多步法需要的初值由经典RK法提供。
f=inline('
-y+2*cos(x)'
'
n=20;
Y=zeros(1,n+1);
Y
(1)=1;
h=0.05*pi;
X=0:
h:
pi;
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,'
gridon;
经典RK法解非线性方程'
legend('
精确解'
经典RK法'
四阶Adams预测-校正算法
h=pi/n;
Y1=Y;
fori=4:
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);
步长四阶Adams预测-校正算法精确解'
Y2'
plot(X,Y2,'
X,Y1,'
四阶Adams预测-校正算法解微分方程'
精确解'
四阶Adams预测-校正算法'