数值分析上机作业总文档格式.docx
- 文档编号:5085777
- 上传时间:2023-05-04
- 格式:DOCX
- 页数:22
- 大小:45.91KB
数值分析上机作业总文档格式.docx
《数值分析上机作业总文档格式.docx》由会员分享,可在线阅读,更多相关《数值分析上机作业总文档格式.docx(22页珍藏版)》请在冰点文库上搜索。
1.0175
0.5073
0.2506
0.1194
0.0477
二、解线性方程组直接法(教材49页15题)
程序如下:
functiontiaojianshu(n)
A=zeros(n);
forj=1:
1:
n
fori=1:
A(i,j)=(1+0.1*i)^(j-1);
c=cond(A)
d=rcond(A)
当n=5时
c=
5.3615e+005
d=
9.4327e-007
当n=10时
8.6823e+011
5.0894e-013
当n=20时
3.4205e+022
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};
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)
Warning:
迭代次数太多,可能不收敛'
本题主程序如下:
functionyakebidiedai
A=[101234;
19-12-3;
2-173-5;
32312-1;
4-3-5-115];
b=[12;
-27;
14;
-17;
12];
x0=[0;
[x,n]=jacobi(A,b,x0)
1.0000
-2.0000
3.0000
n=
67
经过67次迭代,得到最终结果
(2)用Gauss-Seidel迭代法求:
Gauss-Seidel迭代法程序如下:
function[x,n]=gauseidel(A,b,x0,eps,M)
elseifnargin==4
G=(D-L)\U;
f=(D-L)\b;
x=G*x0+f;
x=G*x0+f;
functiongaosidiedai
[x,n]=gauseidel(A,b,x0)
38
经过38次迭代,得到最终结果。
四、矩阵特征值与特征向量的计算(教材100页13题)
幂法求最大特征值的程序:
function[l,v,s]=pmethod(A,x0,eps)
ifnargin==2
eps=1.0e-6;
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)<
eps)
l=m;
s=k;
else
if(k==M)
收敛速度过慢'
s=M;
求解本题程序如下:
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
0
s=
114
结论:
经过114次迭代,求得此矩阵的最大特征值为343.0000,及其对应特征向量为[-0.6667;
-2.0000;
1.0000]
五、函数逼近(教材164页16题)
本题采用最小二乘法进行拟合:
线性最小二乘法程序如下:
function[a,b]=LZXEC(x,y)
if(length(x)==length(y))
n=length(x);
else
x和y的维数不相等'
A=zeros(2,2);
A(2,2)=n;
B=zeros(2,1);
fori=1:
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);
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;
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));
本题计算程序如下(采用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)'
a=a+1;
b=b+1;
plot(q1,q2);
七、非线性方程及非线性方程组的解法
本题采用牛顿法进行求解:
牛顿法解非线性方程程序如下:
function[r,n]=mulNewton(F,x0,eps)
ifnargin==2
eps=1.0e-4;
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;
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);
100000)
迭代步数太多,可能不收敛'
本题解决方案如下:
首先,绘制此方程的图形,大概确定其与X轴的交点位置。
由于
,可以得出
因此绘制程序如下:
ezplot('
log((513+0.6651*x)/(513-0.6651*x))-x/(1400*0.0918)'
[-772,772,-10,10]);
得到图形如下图所示:
经过放大后,发现图形与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=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;
经典RK算法程序如下:
functiony=DELGKT4_lungkuta(f,h,a,b,y0,varvec)
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;
其中FUNVAL函数程序如下:
functionfv=Funval(f,varvec,varval)
iflength(var)<
4
ifvar
(1)==varvec
(1)
fv=subs(f,varvec
(1),varval
(1));
fv=subs(f,varvec
(2),varval
(2));
fv=subs(f,varvec,varval);
本题解决方案如下(步长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.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)