数值分析上机作业总.docx
- 文档编号:2716809
- 上传时间:2023-05-04
- 格式:DOCX
- 页数:22
- 大小:45.91KB
数值分析上机作业总.docx
《数值分析上机作业总.docx》由会员分享,可在线阅读,更多相关《数值分析上机作业总.docx(22页珍藏版)》请在冰点文库上搜索。
数值分析上机作业总
数值分析上机实验
一、解线性方程组直接法(教材49页14题)
追赶法程序如下:
functionx=followup(A,b)
n=rank(A);
for(i=1:
n)
if(A(i,i)==0)
disp('Error:
对角有元素为0');
return;
end
end;
d=ones(n,1);
a=ones(n-1,1);
c=ones(n-1);
for(i=1:
n-1)
a(i,1)=A(i+1,i);
c(i,1)=A(i,i+1);
d(i,1)=A(i,i);
end
d(n,1)=A(n,n);
for(i=2:
n)
d(i,1)=d(i,1)-(a(i-1,1)/d(i-1,1))*c(i-1,1);
b(i,1)=b(i,1)-(a(i-1,1)/d(i-1,1))*b(i-1,1);
end
x(n,1)=b(n,1)/d(n,1);
for(i=(n-1):
-1:
1)
x(i,1)=(b(i,1)-c(i,1)*x(i+1,1))/d(i,1);
end
主程序如下:
functionzhunganfa
A=[2-2000000;-25-200000;0-25-20000;00-25-2000;000-25-200;0000-25-20;00000-25-2;000000-25];
b=[220/27;0;0;0;0;0;0;0];
x=followup(A,b)
计算结果:
x=
8.1478
4.0737
2.0365
1.0175
0.5073
0.2506
0.1194
0.0477
二、解线性方程组直接法(教材49页15题)
程序如下:
functiontiaojianshu(n)
A=zeros(n);
forj=1:
1:
n
fori=1:
1:
n
A(i,j)=(1+0.1*i)^(j-1);
end
end
c=cond(A)
d=rcond(A)
当n=5时
c=
5.3615e+005
d=
9.4327e-007
当n=10时
c=
8.6823e+011
d=
5.0894e-013
当n=20时
c=
3.4205e+022
d=
8.1226e-024
备注:
对于病态矩阵A来说,d为接近0的数;对于非病态矩阵A来说,d为接近1的数。
三、解线性方程组的迭代法(教材74页14题)
(1)用Jacobi迭代法求:
Jacobi迭代法程序如下:
function[x,n]=jacobi(A,b,x0,eps,varargin)
ifnargin==3
eps=1.0e-6;
M=200;
elseifnargin<3
error
return
elseifnargin==5
M=varargin{1};
end
D=diag(diag(A));
L=-tril(A,-1);
U=-triu(A,1);
B=D\(L+U);
f=D\b;
x=B*x0+f;
n=1;
whilenorm(x-x0)>=eps
x0=x;
x=B*x0+f;
n=n+1;
if(n>=M)
disp('Warning:
迭代次数太多,可能不收敛');
return;
end
end
本题主程序如下:
functionyakebidiedai
A=[101234;19-12-3;2-173-5;32312-1;4-3-5-115];
b=[12;-27;14;-17;12];
x0=[0;0;0;0;0];
[x,n]=jacobi(A,b,x0)
计算结果:
x=
1.0000
-2.0000
3.0000
-2.0000
1.0000
n=
67
经过67次迭代,得到最终结果
(2)用Gauss-Seidel迭代法求:
Gauss-Seidel迭代法程序如下:
function[x,n]=gauseidel(A,b,x0,eps,M)
ifnargin==3
eps=1.0e-6;
M=200;
elseifnargin==4
M=200;
elseifnargin<3
error
return;
end
D=diag(diag(A));
L=-tril(A,-1);
U=-triu(A,1);
G=(D-L)\U;
f=(D-L)\b;
x=G*x0+f;
n=1;
whilenorm(x-x0)>=eps
x0=x;
x=G*x0+f;
n=n+1;
if(n>=M)
disp('Warning:
迭代次数太多,可能不收敛');
return;
end
end
本题主程序如下:
functiongaosidiedai
A=[101234;19-12-3;2-173-5;32312-1;4-3-5-115];
b=[12;-27;14;-17;12];
x0=[0;0;0;0;0];
[x,n]=gauseidel(A,b,x0)
计算结果:
x=
1.0000
-2.0000
3.0000
-2.0000
1.0000
n=
38
经过38次迭代,得到最终结果。
四、矩阵特征值与特征向量的计算(教材100页13题)
幂法求最大特征值的程序:
function[l,v,s]=pmethod(A,x0,eps)
ifnargin==2
eps=1.0e-6;
end
v=x0;
M=5000;
m=0;
l=0;
for(k=1:
M)
y=A*v;
m=max(y);
v=y/m;
if(abs(m-l) l=m; s=k; return; else if(k==M) disp('收敛速度过慢'); l=m; s=M; else l=m; end end end 求解本题程序如下: functionmifa A=[19066-8430;6630342-36;336-168147-112;30-3628291]; x0=[0001]'; [l,v,s]=pmethod(A,x0) 求解结果: l= 343.0000 v= -0.6667 -2.0000 0 1.0000 s= 114 结论: 经过114次迭代,求得此矩阵的最大特征值为343.0000,及其对应特征向量为[-0.6667;-2.0000;0;1.0000] 五、函数逼近(教材164页16题) 本题采用最小二乘法进行拟合: 线性最小二乘法程序如下: function[a,b]=LZXEC(x,y) if(length(x)==length(y)) n=length(x); else disp('x和y的维数不相等'); return; end A=zeros(2,2); A(2,2)=n; B=zeros(2,1); fori=1: n A(1,1)=A(1,1)+x(i)*x(i); A(1,2)=A(1,2)+x(i); B(1,1)=B(1,1)+x(i)*y(i); B(2,1)=B(2,1)+y(i); end A(2,1)=A(1,2); s=A\B; a=s (1); b=s (2); 首先利用1/y代替y,1/x代替x并采用线性最小二乘法求出a与b: functionzuixiaoercheng x1=[2356791011121416171920]; y1=[106.42108.26109.58109.50109.86110.00109.93110.59110.60110.72110.90110.76111.10111.30]; x2=1./x1; y2=1./y1; [b,a]=LZXEC(x2,y2) 计算结果: b= 8.4169e-004 a= 0.0090 绘制图形: fplot('x/(0.0090*x+8.4169e-004)',[2,20]) gridon title('最小二乘拟合') 六、数值微分与数值积分(教材207页26题) 本题采用高斯—勒让德求积公式求解: 高斯—勒让德求积公式程序如下: functionI=IntGauss(f,a,b,n,AK,XK) if(n<5&&nargin==4) AK=0; XK=0; else XK1=((b-a)/2)*XK+((a+b)/2); I=((b-a)/2)*sum(AK.*subs(sym(f),findsym(f),XK1)); End ta=(b-a)/2; tb=(a+b)/2; switchn case0, I=2*ta*subs(sym(f),findsym(sym(f)),tb); case1, I=ta*(subs(sym(f),findsym(sym(f)),ta*0.5773503+tb)+... subs(sym(f),findsym(sym(f)),-ta*0.5773503+tb)); case2, I=ta*(0.55555556*subs(sym(f),findsym(sym(f)),ta*0.7745967+tb)+... 0.55555556*subs(sym(f),findsym(sym(f)),-ta*0.7745967+tb)+... 0.88888889*subs(sym(f),findsym(sym(f)),tb)); case3, I=ta*(0.3478548*subs(sym(f),findsym(sym(f)),ta*0.8611363+tb)+... 0.3478548*subs(sym(f),findsym(sym(f)),-ta*0.8611363+tb)+... 0.6521452*subs(sym(f),findsym(sym(f)),ta*0.3398810+tb)... +0.6521452*subs(sym(f),findsym(sym(f)),-ta*0.3398810+tb)); case4, I=ta*(0.2369269*subs(sym(f),findsym(sym(f)),ta*0.9061793+tb)+... 0.2369269*subs(sym(f),findsym(sym(f)),-ta*0.9061793+tb)+... 0.4786287*subs(sym(f),findsym(sym(f)),ta*0.5384693+tb)... +0.4786287*subs(sym(f),findsym(sym(f)),-ta*0.5384693+tb)+... 0.5688889*subs(sym(f),findsym(sym(f)),tb)); end 本题计算程序如下(采用4个节点): functionshuzhijifen a=1; b=1; fors=-5: 0.05: 5 q1(1,a)=IntGauss('cos(1/2*x^2)',0,s,4); q2(1,b)=IntGauss('sin(1/2*x^2)',0,s,4); a=a+1; b=b+1; end plot(q1,q2); 绘制图形: 七、非线性方程及非线性方程组的解法 本题采用牛顿法进行求解: 牛顿法解非线性方程程序如下: function[r,n]=mulNewton(F,x0,eps) ifnargin==2 eps=1.0e-4; end x0=transpose(x0); Fx=subs(F,findsym(F),x0); var=findsym(F); dF=Jacobian(F,var); dFx=subs(dF,findsym(dF),x0); r=x0-inv(dFx)*Fx; n=1; tol=1; whiletol>eps x0=r; Fx=subs(F,findsym(F),x0); dFx=subs(dF,findsym(dF),x0); r=x0-inv(dFx)*Fx; tol=norm(r-x0); n=n+1; if(n>100000) disp('迭代步数太多,可能不收敛'); return; end end 本题解决方案如下: 首先,绘制此方程的图形,大概确定其与X轴的交点位置。 由于 ,可以得出 因此绘制程序如下: ezplot('log((513+0.6651*x)/(513-0.6651*x))-x/(1400*0.0918)',[-772,772,-10,10]); gridon 得到图形如下图所示: 经过放大后,发现图形与x轴的交点接近 处。 计算非零根: 令 为牛顿法接非线性方程的初值。 程序如下: symsx f=log((513+0.6651*x)/(513-0.6651*x))-x/(1400*0.0918); x0=-765 [r,n]=mulNewton(f,x0) 解得: r=-767.3861 x0=765 [r,n]=mulNewton(f,x0) 解得: r=767.3861 结论: 此方程的两个非零根分别为: 八、常微分方程数值解法(教材266页19题) 本题分别采用四阶ADAMS预测—校正算法和经典RK法进行求解: 四阶ADAMS预测—校正算法如下: functiony=DEYCJZ_yds(f,h,a,b,y0,varvec) formatlong; N=(b-a)/h; y=zeros(N+1,1); x=a: h: b; y (1)=y0; y (2)=y0+h*Funval(f,varvec,[x (1)y (1)]); y(3)=y (2)+h*Funval(f,varvec,[x (2)y (2)]); y(4)=y(3)+h*Funval(f,varvec,[x(3)y(3)]); fori=5: N+1 v1=Funval(f,varvec,[x(i-4)y(i-4)]); v2=Funval(f,varvec,[x(i-3)y(i-3)]); v3=Funval(f,varvec,[x(i-2)y(i-2)]); v4=Funval(f,varvec,[x(i-1)y(i-1)]); t=y(i-1)+h*(55*v4-59*v3+37*v2-9*v1)/24; ft=Funval(f,varvec,[x(i)t]); y(i)=y(i-1)+h*(9*ft+19*v4-5*v3+v2)/24; end 经典RK算法程序如下: functiony=DELGKT4_lungkuta(f,h,a,b,y0,varvec) formatlong; N=(b-a)/h; y=zeros(N+1,1); y (1)=y0; x=a: h: b; var=findsym(f); fori=2: N+1 K1=Funval(f,varvec,[x(i-1)y(i-1)]); K2=Funval(f,varvec,[x(i-1)+h/2y(i-1)+K1*h/2]); K3=Funval(f,varvec,[x(i-1)+h/2y(i-1)+K2*h/2]); K4=Funval(f,varvec,[x(i-1)+hy(i-1)+h*K3]); y(i)=y(i-1)+h*(K1+2*K2+2*K3+K4)/6; end 其中FUNVAL函数程序如下: functionfv=Funval(f,varvec,varval) var=findsym(f); iflength(var)<4 ifvar (1)==varvec (1) fv=subs(f,varvec (1),varval (1)); else fv=subs(f,varvec (2),varval (2)); end else fv=subs(f,varvec,varval); end 本题解决方案如下(步长h=0.1): 程序: functionchangweifen symsxy; z=-y+2*cos(x); yy1=DELGKT4_lungkuta(z,0.1,0,pi,1,[xy]); yy2=DEYCJZ_yds(z,0.1,0,pi,1,[xy]); a=0: 0.1: pi; yy3=cos(a)+sin(a); plot(a,yy1,'r*'); gridon; holdon; plot(a,yy2,'b+'); plot(a,yy3,'-go'); legend('RK','Adams','准确解') 整体图形: 局部放大图形: 继续放大: 数据结果: yy1(RK) yy2(ADAMS) yy3(准确解) yy1-yy3 yy2-yy3 1 1 1 0 0 1.094837464 1.1 1.094837582 -1.18389E-07 0.005162418 1.178735678 1.189000833 1.178735909 -2.30293E-07 0.010264924 1.250856361 1.266114065 1.250856696 -3.351E-07 0.01525737 1.310478904 1.324215642 1.310479336 -4.32217E-07 0.013736305 1.357007579 1.369446642 1.3570081 -5.21089E-07 0.012438542 1.389977487 1.401242287 1.389978088 -6.012E-07 0.011264199 1.409059202 1.419252042 1.409059875 -6.72089E-07 0.010192167 1.414062067 1.423285143 1.4140628 -7.33353E-07 0.009222343 1.404936093 1.41328163 1.404936878 -7.84658E-07 0.008344753 1.381772465 1.389323897 1.381773291 -8.25739E-07 0.007550607 1.344802625 1.351635461 1.344803481 -8.56415E-07 0.00683198 1.294395964 1.300578527 1.29439684 -8.76584E-07 0.006181686 1.231056128 1.236650238 1.231057014 -8.86229E-07 0.005593224 1.155415987 1.160477584 1.155416873 -8.85423E-07 0.005060711 1.068231314 1.072811013 1.068232188 -8.74325E-07 0.004578825 0.970373228 0.974516833 0.970374081 -8.53183E-07 0.004142752 0.862819494 0.866568452 0.862820316 -8.22334E-07 0.003748136 0.746644754 0.75003657 0.746645536 -7.82199E-07 0.003391034 0.623009788 0.626078402 0.623010521 -7.33279E-07 0.003067882 0.493149914 0.495926042 0.49315059 -6.76156E-07 0.002775452 0.358362651 0.360874088 0.358363262 -6.11484E-07 0.002510826 0.219994747 0.222266651 0.219995287 -5.39985E-07 0.002271364 0.079428728 0.081483866 0.079429191 -4.62442E-07 0.002054676 -0.061930915 -0.060071936 -0.061930535 -3.7969E-07 0.001858599 -0.202671764 -0.200990293 -0.202671471 -2.92614E-07 0.001681178 -0.341387584 -0.339866738 -0.341387382 -2.02132E-07 0.001520644 -0.476692371 -0.475316868 -0.476692262 -1.09196E-07 0.001375394 -0.607234205 -0.605990211 -0.607234191 -1.47751E-08 0.00124398 -0.731708756 -0.730583746 -0.731708836 8.01498E-08 0.00112509 -0.848872314 -0.847854952 -0.848872489 1.74596E-07 0.001017537 -0.95755422 -0.95663424 -0.957554488 2.6759E-07 0.000920247 结论: 通过图形和数据结果可以看出,在本题中利用经典RK方法获得解的精确度比4阶ADAMS方法要高很多。 《数值计算方法》 上机实验 姓名: 刘天博 学号:
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 数值 分析 上机 作业
![提示](https://static.bingdoc.com/images/bang_tan.gif)