数值分析Matlab作业.docx
《数值分析Matlab作业.docx》由会员分享，可在线阅读，更多相关《数值分析Matlab作业.docx（24页珍藏版）》请在冰豆网上搜索。
数值分析Matlab作业
数值分析编程作业
2012年12月
第二章
14.考虑梯形电阻电路的设计，电路如下：
电路中的各个电流{i1,i2,…,i8}须满足下列线性方程组：
这是一个三对角方程组。
设V=220V，R=27
，运用追赶法，求各段电路的电流量。
Matlab程序如下：
functionchase（）%追赶法求梯形电路中各段的电流量
a=input（'请输入下主对角线向量a='）;
b=input（'请输入主对角线向量b='）;
c=input（'请输入上主对角线向量c='）;
d=input（'请输入右端向量d='）;
n=input（'请输入系数矩阵维数n='）;
u
（1）=b
（1）;
fori=2:
n
l（i）=a（i）/u（i-1）;
u（i）=b（i）-c（i-1）*l（i）;
end
y
（1）=d
（1）;
fori=2:
n
y（i）=d（i）-l（i）*y（i-1）;
end
x（n）=y（n）/u（n）;
i=n-1;
whilei>0
x（i）=（y（i）-c（i）*x（i+1））/u（i）;
i=i-1;
end
x
输入如下:
>>chase
请输入下主对角线向量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.14784.07372.03651.01750.50730.25060.11940.0477
第三章
14.试分别用
（1）Jacobi迭代法；
（2）Gauss-Seidel迭代法解线性方程组
迭代初始向量
。
（1）雅可比迭代法程序如下：
functionjacobi（）%Jacobi迭代法
a=input（'请输入系数矩阵a='）;
b=input（'请输入右端向量b='）;
x0=input（'请输入初始向量x0='）;
n=input（'请输入系数矩阵阶数n='）;
er=input（'请输入允许误差er='）;
N=input（'请输入最大迭代次数N='）;
fori=1:
n
forj=1:
n
ifi==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;
whilek<=N%进行迭代
fori=1:
5
ifmax（abs（x（i）-x0（i）））>er
x=m*x+g;
k=k+1;
else
x
return
end
end
continue
end
x
程序执行如下：
>>jacobi
请输入系数矩阵a=[101234;19-12-3;2-173-5;32312-1;4-3-5-115]
请输入右端向量b=[12-2714-1712]'
请输入初始向量x0=[00000]'
请输入系数矩阵阶数n=5
请输入允许误差er=1.0e-6
请输入最大容许迭代次数N=60
x=
1.0000
-2.0000
3.0000
-2.0000
1.0000
（2）高斯－赛德尔迭代法程序如下：
functiongs_sdl（）%gauss-seiddel迭代法
a=input（'请输入系数矩阵a='）;
b=input（'请输入右端向量b='）;
x0=input（'请输入初始向量x0='）;
n=input（'请输入系数矩阵阶数n='）;
er=input（'请输入允许误差er='）;
N=input（'请输入最大迭代次数N='）;
fori=1:
n
forj=1:
n
ifi<=j
l（i,j）=0;
else
l（i,j）=-a（i,j）;
end
end
end
fori=1:
n
forj=1:
n
ifiu（i,j）=-a（i,j）;
else
u（i,j）=0;
end
end
end
fori=1:
n
forj=1:
n
ifi==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;
whilek<=N
fori=1:
5
ifmax（abs（x（i）-x0（i）））>er
x=m*x+g;
k=k+1;
else
x
return
end
end
continue
end
x
执行结果如下：
>>gs_sdl
请输入系数矩阵a=[101234;19-12-3;2-173-5;32312-1;4-3-5-115]
请输入右端向量b=[12-2714-1712]'
请输入初始向量x0=[00000]'
请输入系数矩阵阶数n=5
请输入允许误差er=1.0e-6
请输入最大容许迭代次数N=60
x=
1.0000
-2.0000
3.0000
-2.0000
1.0000
第四章
已知如下矩阵，试用幂法求按模最大的特征值与特征向量。
Matlab程序代码如下：
functionmifa（）
A=input（'请输入系数矩阵A='）;
x0=input（'请输入初始列向量x0='）;
n=input（'请输入向量维数n='）;
er=input（'请输入允许误差er='）;
N=input（'请输入最大容许迭代次数N='）;
k=1;
mu=0;
whilek<=N
fort=1:
n
ifabs（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=[19066-8430;6630342-36;336-168147-112;30-3628291]
请输入初始列向量x0=[0001]'
请输入向量维数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[ty]=func5（n）
x0=linspace（-5,5,n+1）';
y0=1./（1.+4.*x0.^2）;
b=zeros（1,n+1）;
fori=1:
n+1
s=0;
forj=1:
i
t=1;
fork=1:
i
ifk~=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）;
fori=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次插值运行结果：
[bY]=func5（10）
b=
Columns1through4
-0.00000.00000.0027-0.0000
Columns5through8
-0.0514-0.00000.3920-0.0000
Columns9through11
-1.14330.00001.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=
Columns1through12
2.70033.99944.35154.09743.49262.72371.92111.17150.52740.0154-0.3571-0.5960
Columns13through24
-0.7159-0.7368-0.6810-0.5709-0.4278-0.2704-0.11470.02700.14580.23600.29490.3227
Columns25through36
0.32170.29580.25040.19150.12550.0588-0.0027-0.0537-0.0900-0.1082-0.1062-0.0830
Columns37through48
-0.03900.02450.10520.20000.30500.41580.52800.63690.73790.82690.90020.9549
Columns49through60
0.98861.00000.98860.95490.90020.82690.73790.63690.52800.41580.30500.2000
Columns61through72
0.10520.0245-0.0390-0.0830-0.1062-0.1082-0.0900-0.0537-0.00270.05880.12550.1915
Columns73through84
0.25040.29580.32170.32270.29490.23600.14580.0270-0.1147-0.2704-0.4278-0.5709
Columns85through96
-0.6810-0.7368-0.7159-0.5960-0.35710.01540.52741.17151.92112.72373.49264.0974
Columns97through99
4.35153.99942.7003
绘制原函数和拟合多项式的图形代码：
plot（x,1./（1+4.*x.^2））
holdall
plot（x,y,'r'）
xlabel（'X'）
ylabel（'Y'）
title（'Runge现象'）
gtext（'原函数'）
gtext（'十次牛顿插值多项式'）
绘制结果：
误差计数并绘制误差图：
holdoff
ey=1./（1+4.*x.^2）-y
ey=
Columns1through12
-2.6900-3.9887-4.3403-4.0857-3.4804-2.7109-1.9077-1.1575-0.5128-0.00000.37330.6130
Columns13through24
0.73390.75580.70100.59210.45020.29430.14010.0000-0.1169-0.2051-0.2617-0.2870
Columns25through36
-0.2832-0.2542-0.2053-0.1424-0.0719-0.00000.06740.12540.16960.19710.20620.1962
Columns37through48
0.16790.12340.06600.0000-0.0691-0.1349-0.1902-0.2270-0.2379-0.2171-0.1649-0.0928
Columns49through60
-0.02710-0.0271-0.0928-0.1649-0.2171-0.2379-0.2270-0.1902-0.1349-0.06910.0000
Columns61through72
0.06600.12340.16790.19620.20620.19710.16960.12540.06740.0000-0.0719-0.1424
Columns73through84
-0.2053-0.2542-0.2832-0.2870-0.2617-0.2051-0.11690.00000.14010.29430.45020.5921
Columns85through96
0.70100.75580.73390.61300.37330.0000-0.5128-1.1575-1.9077-2.7109-3.4804-4.0857
Columns97through99
-4.3403-3.9887-2.6900
plot（x,ey）
xlabel（'X'）
ylabel（'ey'）
title（'Runge现象误差图'）
第六章
16、钢包问题。
炼钢唱出钢时所用的盛钢水的钢包，在使用过程中由于钢液及炉渣对包衬耐火材料的侵蚀，使其容积不断增大．经实验，钢包的容积与相应的使用次数的数据如下：
使用次数x容积y
使用次数x容积y
2106.42
3108.26
5109.58
6109.50
7109.86
9110.00
10109.93
11110.59
12110.60
14110.72
16110.90
17110.76
19111.10
20111.30
选用双曲线
对数据进行拟合，使用最小二乘法拟合．
Matlab程序如下：
functiona=nihehanshu（）
x0=[2356791011121416171920];
y0=[106.42108.26109.58109.50109.86110.00109.93110.59110.60110.72110.90110.76111.10111.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'）;
holdon;
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.00900.0008
则拟合曲线为
拟合曲线图、散点图、误差图如下：
第七章
26.考纽螺线的形状像钟表的发条，也称回旋曲线，它在直角坐标系中的参数方程为
曲线关于原点对称。
取a=1,参数s的变化范围[-5，5]，容许误差限分别是10-3和10-7。
选取适当的节点个数，利用数值积分方法计算曲线上点的坐标，并画出曲线的图形。
程序代码如下所示：
functionhuixuan（）%用梯形公式的逐次分半算法计算回旋曲线上点的坐标
er=input（'请选择允许误差1.0e-3或1.0e-7：
'）;
i=1;%x向量分量的下标
fors=-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;
whileabs（T-T0）>3*er
Fx=0;
T0=T;
fork=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向量分量的下标
fors=-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;
whileabs（T-T0）>3*er
Fy=0;
T0=T;
fork=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'）;
ifer==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.求方程
在
附近的根，精确到
。
（1）取
，用简单迭代法
计算；
（2）用加快收敛的迭代格式
，
计算。
程序及计算过程如下：
建一M文件f.m存储函数，
functionf=f（x）
f=exp（-x）;
取
用简单迭代法
计算，Matlab程序如下：
function[x,i]=diedai1（x0）
x=f（x0）;
i=1;
y（i）=x;
whileabs（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次迭代达到精度要求。
用加速收敛法的迭代格式
计算，程序如下：
function[x,i]=diedai2（x0）
w=0.625;
x=w*f（x0）+（1-w）*x0;
i=1;
y（i）=x;
whileabs（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
结果比较
简单迭代法和加速迭代格式的比较
i
x
简单迭代法
30
0.56714328761116
加速的迭代格式
5
0.567143290310401
可见，加速迭代格式收敛比简单迭代格式快。
第九章
设有常微分方程初值问题
其精确解为
。
选取步长使四阶Adams预测-校正算法和经典RK法均稳定，分别用这两种方法求解微分方程，将数值解与精确解进行比较，输出结果。
其中多步法需要的初值由经典RK法提供。
（1）用经典四阶RK法求解，程序代码如下：
functionclassic_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;
fork=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预测-校正算法求解，程序代码如下：
functionadams4（）
n=input（'请输入插值节点数n='）;
h=（pi-0）/n;
x=0:
h:
pi;
fork=1:
n+1
f0（k）=cos（x（k））+sin（x（k））;%f0=cosx+sinx为精确值
end
y
（1）=1;
fork=2:
4%用四阶RK法获得起步值
升级会员 