数值分析Matlab作业Word文档格式.docx
- 文档编号:19775711
- 上传时间:2023-01-10
- 格式:DOCX
- 页数:24
- 大小:246.52KB
数值分析Matlab作业Word文档格式.docx
《数值分析Matlab作业Word文档格式.docx》由会员分享,可在线阅读,更多相关《数值分析Matlab作业Word文档格式.docx(24页珍藏版)》请在冰豆网上搜索。
运行结果如下:
x=8.14784.07372.03651.01750.50730.25060.11940.0477
第三章
14.试分别用
(1)Jacobi迭代法;
(2)Gauss-Seidel迭代法解线性方程组
迭代初始向量
。
(1)雅可比迭代法程序如下:
functionjacobi()%Jacobi迭代法
请输入系数矩阵a='
请输入右端向量b='
x0=input('
请输入初始向量x0='
请输入系数矩阵阶数n='
er=input('
请输入允许误差er='
N=input('
请输入最大迭代次数N='
fori=1:
forj=1:
ifi==j
d(i,j)=a(i,j);
else
d(i,j)=0;
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;
x
return
continue
程序执行如下:
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)高斯-赛德尔迭代法程序如下:
functiongs_sdl()%gauss-seiddel迭代法
ifi<
=j
l(i,j)=0;
l(i,j)=-a(i,j);
j
u(i,j)=-a(i,j);
u(i,j)=0;
m=(d-l)\u;
g=(d-l)\b;
=N
执行结果如下:
gs_sdl
x=
1.0000
第四章
已知如下矩阵,试用幂法求按模最大的特征值与特征向量。
Matlab程序代码如下:
functionmifa()
A=input('
请输入系数矩阵A='
请输入初始列向量x0='
请输入向量维数n='
请输入最大容许迭代次数N='
mu=0;
fort=1:
ifabs(x0(t))==max(abs(x0))
alfa=x0(t);
xb=t;
%最大的x0(i)的下标
y=x0./alfa;
x0=A*y;
lamda=x0(xb);
lamda%按模最大的特征值
x0%按模最大的特征值对应的特征向量
程序执行结果如下:
mifa
请输入系数矩阵A=[19066-8430;
6630342-36;
336-168147-112;
30-3628291]
请输入初始列向量x0=[0001]'
请输入向量维数n=4
请输入最大容许迭代次数N=100
lamda=343.0000
x0=
114.3333
343.0000
-0.0000
-171.5002
第五章
试编写MATLAB函数实现Newton插值,要求能输出插值多项式。
对函数
在区间[-5,5]上实现10次多项式插值。
%此函数实现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);
n+1
s=0;
i
t=1;
fork=1:
ifk~=j
t=(x0(j)-x0(k))*t;
end;
s=s+y0(j)/t;
b(i)=s;
end;
t=linspace(0,0,n+1);
s=linspace(0,0,n+1);
s(n+1-i:
n+1)=b(i+1).*poly(x0(1:
i));
t=t+s;
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('
原函数'
十次牛顿插值多项式'
绘制结果:
误差计数并绘制误差图:
holdoff
ey=1./(1+4.*x.^2)-y
ey=
-2.6900-3.9887-4.3403-4.0857-3.4804-2.7109-1.9077-1.1575-0.5128-0.00000.37330.6130
0.73390.75580.70100.59210.45020.29430.14010.0000-0.1169-0.2051-0.2617-0.2870
-0.2832-0.2542-0.2053-0.1424-0.0719-0.00000.06740.12540.16960.19710.20620.1962
0.16790.12340.06600.0000-0.0691-0.1349-0.1902-0.2270-0.2379-0.2171-0.1649-0.0928
-0.02710-0.0271-0.0928-0.1649-0.2171-0.2379-0.2270-0.1902-0.1349-0.06910.0000
0.06600.12340.16790.19620.20620.19710.16960.12540.06740.0000-0.0719-0.1424
-0.2053-0.2542-0.2832-0.2870-0.2617-0.2051-0.11690.00000.14010.29430.45020.5921
0.70100.75580.73390.61300.37330.0000-0.5128-1.1575-1.9077-2.7109-3.4804-4.0857
-4.3403-3.9887-2.6900
plot(x,ey)
xlabel('
ylabel('
ey'
title('
Runge现象误差图'
第六章
16、钢包问题。
炼钢唱出钢时所用的盛钢水的钢包,在使用过程中由于钢液及炉渣对包衬耐火材料的侵蚀,使其容积不断增大.经实验,钢包的容积与相应的使用次数的数据如下:
使用次数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
选用双曲线
对数据进行拟合,使用最小二乘法拟合.
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-'
拟合曲线误差'
subplot(1,2,1);
plot(x0,y0,'
go'
holdon;
x=2:
0.5:
20;
y=1./(a
(1)+a
(2)*1./x);
r*-'
legend('
散点'
'
拟合曲线图1/y=a
(1)+a
(2)*1/x'
最小二乘法拟合曲线'
求的系数为:
0.00900.0008
则拟合曲线为
拟合曲线图、散点图、误差图如下:
第七章
26.考纽螺线的形状像钟表的发条,也称回旋曲线,它在直角坐标系中的参数方程为
曲线关于原点对称。
取a=1,参数s的变化范围[-5,5],容许误差限分别是10-3和10-7。
选取适当的节点个数,利用数值积分方法计算曲线上点的坐标,并画出曲线的图形。
程序代码如下所示:
functionhuixuan()%用梯形公式的逐次分半算法计算回旋曲线上点的坐标
请选择允许误差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;
2^(m-1)%计算新增加节点处的函数值之和
fx3=cos((a+(2*k-1)*h)^2/2);
Fx=Fx+fx3;
T=T0/2+h*Fx;
m=m+1;
h=h/2;
x(i)=T;
i=i+1;
j=1;
%y向量分量的下标
n=1;
fy1=sin(a^2/2);
fy2=sin(b^2/2);
T=h*(fy1+fy2);
Fy=0;
2^(n-1)
fy3=sin((a+(2*k-1)*h)^2/2);
Fy=Fy+fy3;
T=T0/2+h*Fy;
n=n+1;
y(j)=T;
j=j+1;
k*'
x,y,'
k'
ifer==1.0e-3
er=1.0e-3'
else
er=1.0e-7'
huixuan
1.0e-3
1.0e-7
第八章
20.求方程
在
附近的根,精确到
(1)取
,用简单迭代法
计算;
(2)用加快收敛的迭代格式
,
计算。
程序及计算过程如下:
建一M文件f.m存储函数,
functionf=f(x)
f=exp(-x);
取
用简单迭代法
计算,Matlab程序如下:
function[x,i]=diedai1(x0)
x=f(x0);
y(i)=x;
whileabs(x-x0)>
10^-8
x0=x;
x=f(x);
y(i)=x;
取初始值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;
x=w*f(x)+(1-w)*x;
同样取x0=0.5,得
x=
0.567143290310401
i=
结果比较
简单迭代法和加速迭代格式的比较
简单迭代法
0.56714328761116
加速的迭代格式
0.567143290310401
可见,加速迭代格式收敛比简单迭代格式快。
第九章
设有常微分方程初值问题
其精确解为
选取步长使四阶Adams预测-校正算法和经典RK法均稳定,分别用这两种方法求解微分方程,将数值解与精确解进行比较,输出结果。
其中多步法需要的初值由经典RK法提供。
(1)用经典四阶RK法求解,程序代码如下:
functionclassic_rk4()
请输入插值节点数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:
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);
subplot(3,1,1);
plot(x,f0,'
y=cosx+sinx'
subplot(3,1,2);
经典四阶RK法'
subplot(3,1,3);
T=y-f0;
%计算经典四阶RK法的误差
plot(x,T,'
经典四阶RK法的误差'
classic_rk4
请输入插值节点数n=3000
(2)用四阶Adams预测-校正算法求解,程序代码如下:
functionadams4()
h=(pi-0)/n;
f0(k)=cos(x(k))+sin(x(k));
fork=2:
4%用四阶RK法获得起步值
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 数值 分析 Matlab 作业