数值分析Matlab作业.docx
- 文档编号:6757250
- 上传时间:2023-01-10
- 格式:DOCX
- 页数:24
- 大小:246.52KB
数值分析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
ifi u(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法获得起步值
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 数值 分析 Matlab 作业