数值分析Matlab作业

数值分析编程作业
2012年12月
第二章
14.考虑梯形电阻电路的设计，电路如下：
电路中的各个电流{i1,i2,…,i8}须满足下列线性方程组：
12123
234
345
456
567
6787822/25202520252025202520
2520
250
i i V R i i i i i i i i i i i i i i i i i i i i -=-+-=-+-=-+-=-+-=-+-=-+-=-+=
这是一个三对角方程组。

设V=220V ，R=27Ω，运用追赶法，求各段电路的电流量。

Matlab 程序如下：
function chase () %追赶法求梯形电路中各段的电流量 a=input('请输入下主对角线向量a='); b=input('请输入主对角线向量b='); c=input('请输入上主对角线向量c='); d=input('请输入右端向量d='); n=input('请输入系数矩阵维数n='); u(1)=b(1); for i=2:n
l(i)=a(i)/u(i-1); u(i)=b(i)-c(i-1)*l(i); end
y(1)=d(1); for i=2:n
y(i)=d(i)-l(i)*y(i-1); end
x(n)=y(n)/u(n); i=n-1; while i>0
x(i)=(y(i)-c(i)*x(i+1))/u(i); i=i-1; end x
输入如下:
请输入下主对角线向量a=[0,-2,-2,-2,-2,-2,-2,-2]; 请输入主对角线向量b=[2,5,5,5,5,5,5,5];
请输入上主对角线向量c=[-2,-2,-2,-2,-2,-2,-2,0]; 请输入方程组右端向量d=[220/27,0,0,0,0,0,0,0]; 请输入系数矩阵阶数n=8 运行结果如下：
x = 8.1478 4.0737 2.0365 1.0175 0.5073 0.2506 0.1194 0.0477
第三章
14.试分别用(1)Jacobi 迭代法；(2)Gauss-Seidel 迭代法解线性方程组
1234510123412191232721735143231211743511512x x x x x ⎡⎤⎡⎤⎡⎤
⎢⎥⎢⎥⎢⎥---⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥=--⎢⎥⎢⎥⎢⎥
--⎢⎥⎢⎥⎢⎥
⎢⎥⎢⎥⎢⎥---⎣⎦⎣⎦⎣⎦ 迭代初始向量
(0)(0,0,0,0,0)T x =。

（1）雅可比迭代法程序如下：
function jacobi() %Jacobi 迭代法 a=input('请输入系数矩阵a='); b=input('请输入右端向量b='); x0=input('请输入初始向量x0='); n=input('请输入系数矩阵阶数n='); er=input('请输入允许误差er='); N=input('请输入最大迭代次数N='); for i=1:n for j=1:n if i==j
d(i,j)=a(i,j); else
d(i,j)=0; end end end
m=eye(5)-d\a; %迭代矩阵 g=d\b;
x=m*x0+g; k=1;
while k<=N %进行迭代 for i=1:5
if max(abs(x(i)-x0(i))) >er x=m*x+g; k=k+1;
x
return
end
end
continue
end
x
程序执行如下：
>>jacobi
请输入系数矩阵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]'
请输入初始向量x0=[0 0 0 0 0]'
请输入系数矩阵阶数n=5
请输入允许误差er=1.0e-6
请输入最大容许迭代次数N=60
x =
1.0000
-2.0000
3.0000
-2.0000
1.0000
（2）高斯－赛德尔迭代法程序如下：
function gs_sdl() %gauss-seiddel迭代法
a=input('请输入系数矩阵a=');
b=input('请输入右端向量b=');
x0=input('请输入初始向量x0=');
n=input('请输入系数矩阵阶数n=');
er=input('请输入允许误差er=');
N=input('请输入最大迭代次数N=');
for i=1:n
for j=1:n
if i<=j
l(i,j)=0;
else
l(i,j)=-a(i,j);
end
end
end
for i=1:n
for j=1:n
if i<j
u(i,j)=-a(i,j);
else
u(i,j)=0;
end
end
end
for i=1:n
for j=1:n
if i==j
d(i,j)=a(i,j);
else
d(i,j)=0;
end
end
end
m=(d-l)\u; %迭代矩阵
g=(d-l)\b;
x=m*x0+g;
k=1;
while k<=N
for i=1:5
if max(abs(x(i)-x0(i))) >er
x=m*x+g;
k=k+1;
else
x
return
end
end
continue
end
x
执行结果如下：
>> gs_sdl
请输入系数矩阵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]'
请输入初始向量x0=[0 0 0 0 0]'
请输入系数矩阵阶数n=5
请输入允许误差er=1.0e-6
请输入最大容许迭代次数N=60
x =
1.0000
-2.0000
3.0000
-2.0000
1.0000
已知如下矩阵，试用幂法求按模最大的特征值与特征向量。

190668430663034236336168147112303628291-⎡⎤⎢⎥-⎢⎥⎢⎥--⎢⎥-⎣⎦
Matlab 程序代码如下：
function mifa ()
A=input('请输入系数矩阵A='); x0=input('请输入初始列向量x0='); n=input('请输入向量维数n='); er=input('请输入允许误差er=');
N=input('请输入最大容许迭代次数N='); k=1; mu=0;
while k<=N for t=1:n
if abs(x0(t))==max(abs(x0)) alfa=x0(t);
xb=t; %最大的x0（i ）的下标 end end
y=x0./alfa; x0=A*y;
lamda=x0(xb); k=k+1; end
lamda %按模最大的特征值
x0 %按模最大的特征值对应的特征向量 程序执行结果如下： >> mifa
请输入系数矩阵A=[190 66 -84 30;66 303 42 -36;336 -168 147 -112;30 -36 28 291] 请输入初始列向量x0=[0 0 0 1]' 请输入向量维数n=4
请输入允许误差er=1.0e-6
请输入最大容许迭代次数N=100 lamda = 343.0000 x0 =
114.3333 343.0000 -0.0000 -171.5002
试编写MATLAB函数实现Newton插值，要求能输出插值多项式。

对函数在区间[-5，
5]上实现10次多项式插值。

Matlab程序代码如下：
%此函数实现y=1/(1+4*x^2)的n次Newton插值，n由调用函数时指定
%函数输出为插值结果的系数向量（行向量）和插值多项式
function [t y]=func5(n)
x0=linspace(-5,5,n+1)';
y0=1./(1.+4.*x0.^2);
b=zeros(1,n+1);
for i=1:n+1
s=0;
for j=1:i
t=1;
for k=1:i
if k~=j
t=(x0(j)-x0(k))*t;
end;
end;
s=s+y0(j)/t;
end;
b(i)=s;
end;
t=linspace(0,0,n+1);
for i=1:n
s=linspace(0,0,n+1);
s(n+1-i:n+1)=b(i+1).*poly(x0(1:i));
t=t+s;
end;
t(n+1)=t(n+1)+b(1);
y=poly2sym(t);
10次插值运行结果：
[b Y]=func5(10)
b =
Columns 1 through 4
-0.0000 0.0000 0.0027 -0.0000
Columns 5 through 8
-0.0514 -0.0000 0.3920 -0.0000
Columns 9 through 11
-1.1433 0.0000 1.0000
Y =
- (7319042784910035*x^10)/147573952589676412928 + x^9/18446744073709551616 + (256*x^8)/93425 - x^7/1152921504606846976 - (28947735013693*x^6)/562949953421312 - (3*x^5)/72057594037927936 + (36624*x^4)/93425 - (5*x^3)/36028797018963968 - (5148893614132311*x^2)/4503599627370496 + (7*x)/36028797018963968 + 1
b为插值多项式系数向量，Y为插值多项式。

插值近似值：
x1=linspace(-5,5,101);
x=x1(2:100);
y=polyval(b,x)
y =
Columns 1 through 12
2.7003
3.9994
4.3515 4.0974 3.4926 2.7237 1.9211 1.1715 0.5274 0.0154 -0.3571 -0.5960
Columns 13 through 24
-0.7159 -0.7368 -0.6810 -0.5709 -0.4278 -0.2704 -0.1147 0.0270 0.1458 0.2360 0.2949 0.3227
Columns 25 through 36
0.3217 0.2958 0.2504 0.1915 0.1255 0.0588 -0.0027 -0.0537 -0.0900 -0.1082 -0.1062 -0.0830
Columns 37 through 48
-0.0390 0.0245 0.1052 0.2000 0.3050 0.4158 0.5280 0.6369 0.7379 0.8269 0.9002 0.9549
Columns 49 through 60
0.9886 1.0000 0.9886 0.9549 0.9002 0.8269 0.7379 0.6369 0.5280 0.4158 0.3050 0.2000
Columns 61 through 72
0.1052 0.0245 -0.0390 -0.0830 -0.1062 -0.1082 -0.0900 -0.0537 -0.0027 0.0588 0.1255 0.1915
Columns 73 through 84
0.2504 0.2958 0.3217 0.3227 0.2949 0.2360 0.1458 0.0270 -0.1147 -0.2704 -0.4278 -0.5709
Columns 85 through 96
-0.6810 -0.7368 -0.7159 -0.5960 -0.3571 0.0154 0.5274 1.1715 1.9211 2.7237 3.4926 4.0974
Columns 97 through 99
4.3515 3.9994 2.7003
绘制原函数和拟合多项式的图形代码：
plot(x,1./(1+4.*x.^2))
hold all
plot(x,y,'r')
xlabel('X')
ylabel('Y')
title('Runge现象')
gtext('原函数')
gtext('十次牛顿插值多项式')
绘制结果：
误差计数并绘制误差图：
hold off
ey=1./(1+4.*x.^2)-y
ey =
Columns 1 through 12
-2.6900 -3.9887 -4.3403 -4.0857 -3.4804 -2.7109 -1.9077 -1.1575 -0.5128 -0.0000 0.3733 0.6130
Columns 13 through 24
0.7339 0.7558 0.7010 0.5921 0.4502 0.2943 0.1401 0.0000
-0.1169 -0.2051 -0.2617 -0.2870
Columns 25 through 36
-0.2832 -0.2542 -0.2053 -0.1424 -0.0719 -0.0000 0.0674 0.1254
0.1696 0.1971 0.2062 0.1962
Columns 37 through 48
0.1679 0.1234 0.0660 0.0000 -0.0691 -0.1349 -0.1902 -0.2270
-0.2379 -0.2171 -0.1649 -0.0928
Columns 49 through 60
-0.0271 0 -0.0271 -0.0928 -0.1649 -0.2171 -0.2379 -0.2270
-0.1902 -0.1349 -0.0691 0.0000
Columns 61 through 72
0.0660 0.1234 0.1679 0.1962 0.2062 0.1971 0.1696 0.1254
0.0674 0.0000 -0.0719 -0.1424
Columns 73 through 84
-0.2053 -0.2542 -0.2832 -0.2870 -0.2617 -0.2051 -0.1169 0.0000
0.1401 0.2943 0.4502 0.5921 Columns 85 through 96
0.7010 0.7558 0.7339 0.6130 0.3733 0.0000 -0.5128 -1.1575 -1.9077 -2.7109 -3.4804 -4.0857 Columns 97 through 99
-4.3403 -3.9887 -2.6900
plot(x,ey) xlabel('X') ylabel('ey')
title('Runge 现象误差图
')
第六章
16、钢包问题。

炼钢唱出钢时所用的盛钢水的钢包，在使用过程中由于钢液及炉渣对包衬耐火
选用双曲线11
*
a b
y x
=+
对数据进行拟合，使用最小二乘法拟合．
Matlab程序如下：
function a=nihehanshu()
x0=[2 3 5 6 7 9 10 11 12 14 16 17 19 20];
y0=[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];
A=zeros(2,2);
B=zeros(2,1);
a=zeros(2,1);
x=1./x0;
y=1./y0;
A(1,1)=14;
A(1,2)=sum(x);
A(2,1)=A(1,2);
A(2,2)=sum(x.^2);
B(1)=sum(y);
B(2)=sum(x.*y);
a=A\B;
y=1./(a(1)+a(2)*1./x0);
subplot(1,2,2);
plot(x0,y0-y,'bd-');
title('拟合曲线误差');
subplot(1,2,1);
plot(x0,y0,'go');
hold on;
x=2:0.5:20;
y=1./(a(1)+a(2)*1./x);
plot(x,y,'r*-');
legend('散点' ,'拟合曲线图1/y=a(1)+a(2)*1/x');
title('最小二乘法拟合曲线');
求的系数为：0.0090 0.0008
则拟合曲线为
x y
1
0008
.0
009
.0
1
+
=
拟合曲线图、散点图、误差图如下：
第七章
26. 考纽螺线的形状像钟表的发条，也称回旋曲线，它在直角坐标系中的参数方程为
()()⎪⎩⎪⎨⎧
==⎰⎰s s dt
a s y dt a s x t t 020221sin 21cos
曲线关于原点对称。

取a=1,参数s 的变化范围[-5，5]，容许误差限分别是10-3和10-7。

选取适当的节点个数，利用数值积分方法计算曲线上点的坐标，并画出曲线的图形。

程序代码如下所示：
function huixuan () %用梯形公式的逐次分半算法计算回旋曲线上点的坐标 er=input('请选择允许误差1.0e-3或1.0e-7：');
i=1; % x 向量分量的下标 for s=-5:0.1:5 m=1; b=s; a=0;
h=(b-a)/2;
fx1=cos(a^2/2); fx2=cos(b^2/2); T=h*(fx1+fx2); T0=5;
while abs(T-T0)>3*er Fx=0; T0=T;
for k=1:2^(m-1) %计算新增加节点处的函数值之和 fx3=cos((a+(2*k-1)*h)^2/2);
Fx=Fx+fx3;
end
T=T0/2+h*Fx;
m=m+1;
h=h/2;
end
x(i)=T;
i=i+1;
end
j=1; %y向量分量的下标for s=-5:0.1:5
n=1;
b=s;
a=0;
h=(b-a)/2;
fy1=sin(a^2/2);
fy2=sin(b^2/2);
T=h*(fy1+fy2);
T0=5;
while abs(T-T0)>3*er
Fy=0;
T0=T;
for k=1:2^(n-1)
fy3=sin((a+(2*k-1)*h)^2/2);
Fy=Fy+fy3;
end
T=T0/2+h*Fy;
n=n+1;
h=h/2;
end
y(j)=T;
j=j+1;
end
plot(x,y,'k*',x,y,'k');
if er==1.0e-3
title('er=1.0e-3');
else
title('er=1.0e-7');
end
程序执行结果如下：
>> huixuan
请选择允许误差1.0e-3或1.0e-7：1.0e-3
>> huixuan
请选择允许误差1.0e-3或1.0e-7：1.0e-7
第八章
20.求方程
x
x e-
=在0.5
x=附近的根，精确到8
10-。

（1）取
()x
x e
ϕ-
=
，用简单迭代法1
()
n n
x x
ϕ
+
=
计算；
（2） 用加快收敛的迭代格式1()(1)n n n
x x x λϕλ+=+-，0.625λ=计算。

程序及计算过程如下： 建一M 文件f.m 存储函数， function f=f(x) f=exp(-x); 取x
e
x -=)(ϕ,用简单迭代法)(1n n x x ϕ=+计算，Matlab 程序如下：
function [x,i]=diedai1(x0) x=f(x0); i=1; y(i)=x;
while abs(x-x0)>10^-8 i=i+1; x0=x; x=f(x); y(i)=x; end
取初始值x0=0.5,输入[x,i]=diedai1(0.5)得结果 x =
0.567143287611168 i =
30 可以看出用简单收敛法经过30次迭代达到精度要求。

用加速收敛法的迭代格式625.0,)1()(1=-+=+λλλϕn n n x x x 计算，程序如下： function [x,i]=diedai2(x0) w=0.625;
x=w*f(x0)+(1-w)*x0; i=1; y(i)=x;
while abs(x-x0)>10^-8 i=i+1; x0=x;
x=w*f(x)+(1-w)*x; y(i)=x; end
同样取x0=0.5，得 x =
0.567143290310401 i =
5 结果比较
简单迭代法和加速迭代格式的比较
可见，加速迭代格式收敛比简单迭代格式快。

第九章
设有常微分方程初值问题⎩⎨
⎧=<<+-=1)0(0cos 2'y x x y y π
其精确解为x x y sin cos +=。

选取步长使四阶Adams 预测-校正算法和经典RK 法均稳定，分别用这两种方法求解微分方程，将数值解与精确解进行比较，输出结果。

其中多步法需要的初值
由经典RK 法提供。

（1）用经典四阶RK 法求解，程序代码如下： function classic_rk4()
n=input('请输入插值节点数n='); y(1)=1;
f0(1)=1; %f0=cosx+sinx 为精确值 h=pi/n; %步长 x=0:h:pi; k=2;
eps=1.0e-3; for k=1:n
f0(k+1)=cos(x(k))+sin(x(k)); k1=-y(k)+2*cos(x(k));
k2=-(y(k)+h*k1/2)+2*cos(x(k)+h/2); k3=-(y(k)+h*k2/2)+2*cos(x(k)+h/2); k4=-(y(k)+h*k3)+2*cos(x(k)+h);
y(k+1)=y(k)+h/6*(k1+2*k2+2*k3+k4); end
subplot(3,1,1); plot(x,f0,'k');
title('y=cosx+sinx'); subplot(3,1,2); plot(x,y,'k');
title('经典四阶RK 法'); subplot(3,1,3);
T=y-f0; %计算经典四阶RK 法的误差 plot(x,T,'k');
title('经典四阶RK 法的误差'); 程序执行结果如下： >> classic_rk4
请输入插值节点数n=3000
（2）用四阶Adams预测-校正算法求解，程序代码如下：
function adams4()
n=input('请输入插值节点数n=');
h=(pi-0)/n;
x=0:h:pi;
for k=1:n+1
f0(k)=cos(x(k))+sin(x(k)); %f0=cosx+sinx为精确值
end
y(1)=1;
for k=2:4 %用四阶RK法获得起步值
k1=h*(-y(k-1)+2*cos(x(k-1)));
k2=h*(-(y(k-1)+k1/2)+2*cos(x(k-1)+h/2));
k3=h*(-(y(k-1)+k2/2)+2*cos(x(k-1)+h/2));
k4=h*(-(y(k-1)+k3)+2*cos(x(k-1)+h));
y(k)=y(k-1)+h/6*(k1+2*k2+2*k3+k4);
end
f=-y(k)+2*cos(x(k)); %fn
f1=-y(k-1)+2*cos(x(k-1)); %fn-1
f2=-y(k-2)+2*cos(x(k-2)); %fn-2
f3=-y(k-3)+2*cos(x(k-3)); %fn-3
pre=0;
mod=0;
for j=k+1:n+1
pre1=y(j-1)+h/24*(55*f-59*f1+37*f2-9*f3); %计算预测值
l=pre1+251/270*(mod-pre); %用局部截断误差进一步修正预测值mod1=y(j-1)+h/24*(9*(-l+2*cos(x(j)))+19*f-5*f1+f2); %计算校正值
y(j)=mod1-19/270*(mod1-pre1); %用局部截断误差进一步修正校正值f=f1;
f1=f2;
f2=f3;
f3=-y(j)+2*cos(x(j));
pre=pre1;
mod=mod1;
end
subplot(3,1,1);
plot(x,f0,'k');
title('y=cosx+sinx');
subplot(3,1,2);
plot(x,y,'k');
title('Adams预测校正公式法');
T=y-f0; %Adams预测校正公式的误差
subplot(3,1,3);
plot(x,T,'k');
title('Adams预测校正公式的误差');
程序执行结果如下：
>> adams4
请输入插值节点数n=3000。