数值分析课程设计大作业Word格式.docx
- 文档编号:3039189
- 上传时间:2023-05-01
- 格式:DOCX
- 页数:76
- 大小:533.30KB
数值分析课程设计大作业Word格式.docx
《数值分析课程设计大作业Word格式.docx》由会员分享,可在线阅读,更多相关《数值分析课程设计大作业Word格式.docx(76页珍藏版)》请在冰点文库上搜索。
ifp==fix(p)
break
end
disp([x,p])
执行代码后得:
n=102315621(输入n=1000000)
即最后每个人分得1023个椰子,椰子总数为15621
1.2当
时,选择稳定的算法计算积分
由
得
由上式可知求
时,
的误差的影响被缩小了。
n=100时
的近似值为0。
matlab代码为
fprintf('
稳定算法:
\n'
)
y0=0;
n=100;
plot(n,y0,'
r*'
holdon
y[100]=%10.6f'
y0);
while
(1)
y1=1/10*[(1-exp(-n))/n-y0];
fprintf('
y[%10.0f]=%10.6f'
n-1,y1);
plot(n-1,y1,'
)
if(n<
=1)break;
y0=y1;
n=n-1;
ifmod(n,3)==0,fprintf('
),end,end
(具体值已省略)
编程实现得下图。
由图可知,该算法是稳定的。
1.3绘制静态和动态的Koch分形曲线
Koch曲线程序koch.m
functionkoch(a1,b1,a2,b2,n)
%koch(0,0,9,0,3)
%a1,b1,a2,b2为初始线段两端点坐标,n为迭代次数
%例如a1=0;
b1=0;
a2=9;
b2=0;
n=3;
%第i-1次迭代时由各条线段产生的新四条线段的五点横、纵坐标存储在数组A、B中
[A,B]=sub_koch1(a1,b1,a2,b2);
fori=1:
n
forj=1:
length(A)/5;
w=sub_koch2(A(1+5*(j-1):
5*j),B(1+5*(j-1):
5*j));
4
[AA(5*4*(j-1)+5*(k-1)+1:
5*4*(j-1)+5*(k-1)+5),BB(5*4*(j-1)+5*(k-1)+1:
5*4*(j-1)+5*(k-1)+5)]=sub_koch1(w(k,1),w(k,2),w(k,3),w(k,4));
end
end
A=AA;
B=BB;
plot(A,B)
axisequal
%由以(ax,ay),(bx,by)为端点的线段生成新的中间三点坐标并把这五点横、纵坐标依次分别存%储在数组A,B中
function[A,B]=sub_koch1(ax,ay,bx,by)
cx=ax+(bx-ax)/3;
cy=ay+(by-ay)/3;
ex=bx-(bx-ax)/3;
ey=by-(by-ay)/3;
L=sqrt((ex-cx).^2+(ey-cy).^2);
alpha=atan((ey-cy)./(ex-cx));
if(ex-cx)<
alpha=alpha+pi;
dx=cx+cos(alpha+pi/3)*L;
dy=cy+sin(alpha+pi/3)*L;
A=[ax,cx,dx,ex,bx];
B=[ay,cy,dy,ey,by];
%把由函数sub_koch1生成的五点横、纵坐标A,B顺次划分为四组,分别对应四条折线段中
%每条线段两端点的坐标,并依次分别存储在4*4阶矩阵k中,k中第i(i=1,2,3,4)行数字代表第%i条线段两端点的坐标
functionw=sub_koch2(A,B)
a11=A
(1);
b11=B
(1);
a12=A
(2);
b12=B
(2);
a21=A
(2);
b21=B
(2);
a22=A(3);
b22=B(3);
a31=A(3);
b31=B(3);
a32=A(4);
b32=B(4);
a41=A(4);
b41=B(4);
a42=A(5);
b42=B(5);
w=[a11,b11,a12,b12;
a21,b21,a22,b22;
a31,b31,a32,b32;
a41,b41,a42,b42];
%调用函数得到下图
n=5;
i=0;
whilei<
figure(i+1);
koch(0,0,3,0,i)
i=i+1;
pause
(1)
end
将pause
(1)去掉可得静态图
2.1小行星轨道问题
为了确定方程中的5个待定系数,需要将上述5个点的坐标代入上面的方程
得:
将这一包含5个未知数的线性方程组,写成矩阵的形式
=
(1)求解这一线性方程组,即可得到曲线方程的系数
X0=[5360558460628596666268894];
Y0=[606211179169542349268894];
A=zeros(5);
X0
(1);
A(i,1)=X0(i)^2;
A(i,2)=2*X0(i)*Y0(i);
A(i,3)=Y0(i)^2;
A(i,4)=2*X0(i);
A(i,5)=2*Y0(i);
end;
formatlongg;
A
A=
28734960256499070203674784410721012124
3417571600130704868012497004111692022353951253881213142297228743811612571833908
4443822244313204740855187406413332446984
474638323694927664724746383236137788137788
B=[-1-1-1-1-1]'
;
x=A\B
x=
-8.06820280371841e-011
-7.63620099622306e-011
-3.0801152978055e-010
-8.89025159419867e-006
2.02829368401655e-005
(2)用Lu分解法解可得
formatlongg;
A=[28734960256499070203674784410721012124;
3417571600130704868012497004111692022358;
3951253881213142297228743811612571833908;
4443822244313204740855187406413332446984;
474638323694927664724746383236137788137788];
[L,U,flag]=LU_Decom(A),formatlongg;
x=U\(L\B)
flag=
OK
-8.06820280370254e-011
-7.63620099622621e-011
-3.08011529780421e-010
-8.89025159420231e-006
2.02829368401615e-005
jacobi迭代法:
jacobic(A)
因为谱半径不小于1,所以Jacobi迭代序列发散,谱半径SRH和B的所有特征值H如下:
SRH=
4.41963931714337
ans=
-4.41963931714337
1.5453292801696
0.757138763648732
1.11838895650533
0.9987823168197
GSC(A)
因为谱半径不小于1,所以Gauss-Seidel迭代序列发散,谱半径SRH和B的所有特征值H如下:
1.12218280703645
0
0.0914047045360394
1.10703104191813+0.183783907450762i
1.10703104191813-0.183783907450762i
0.99748223605848
2.2
(1)用Gauss列主元消去法、Gauss按比例列主元消去法、Cholesky分解求解下列线性方程组,并彼此互相验证。
(2)判断用Jacobi迭代法、Gauss-Seidel迭代法、SOR法(分别取
)解下列线性方程组的收敛性.若收敛,再用Jacobi迭代法、Gauss-Seidel迭代法、SOR法(分别取
)分别解线性方程组,并比较各种方法的收敛速度.
(3)用Cholesky分解求解下列线性方程组
(1)
A=[1-121;
-130-3;
209-6;
1-3-619];
b=[1357]'
[RA,RB,n,x]=liezy(A,b),[RA,RB,n,x]=bilizy(A,b),cholesky(A,b)
列主元
因为RA=RB=n,所以此方程组有唯一解.
RA=
4
RB=
n=
-8.000000000000005
0.333333333333332
3.666666666666668
2.000000000000000
比例主元
-8.000000000000000
0.333333333333333
3.666666666666667
cholesky分解
S1=
0
003.6666666666666732.000000000000003
00.3333333333333293.6666666666666732.000000000000003
-8.0000000000000200.3333333333333293.6666666666666732.000000000000003
-8.0000000000000200.3333333333333293.6666666666666732.000000000000003
对比可知,三种方法可相互验证。
(2)
x0=ones(4,1);
[]
jacobic(A),GSC(A)
因为谱半径小于1,所以Jcobi迭代序列收敛,谱半径SRH和B的所有特征值H如下:
0.966881024532242
0.536020655954138
-0.903317605697383
-0.599584074788997
因为谱半径小于1,所以Gauss-Seidel迭代序列收敛,谱半径SRH和B的所有特征值H如下:
0.934888367800313
0.281485901205534
-2.38957160758552e-017
b=[1357]'
x0=ones(4,1);
[x,n]=jacobi(A,b,x0)
-7.99997357113847
0.333339142428815
3.66665839089021
1.99999680708368
376
gauseidel法
[x,n]=gauseidel(A,b,x0)
-7.99998673303378
0.333336112109258
3.66666262275452
1.99999846346783
197
SOR法:
w1=0.8;
w2=1.2;
w3=1.3;
w4=1.6;
SORC(A,w1),SORC(A,w2),SORC(A,w3),SORC(A,w4)
w=
0.8
因为谱半径小于1,所以SOR迭代序列收敛,谱半径SRH和B的所有特征值H如下:
0.956498743812875
0.503328665301619
0.050190208900126
.0662********
1.2
0.901948033909749
0.901948033909749
0.0392206356524567
0.00773145469258105+0.212532205268192i
0.00773145469258105-0.212532205268192i
1.3
0.877539492310148
0.877539492310148
0.0946172161251248
-0.0537947284866427+0.307669991725664i
-0.0537947284866427-0.307669991725664i
1.6
0.608878775828679
0.590094939083101+0.0369506212620504i
0.590094939083101-0.0369506212620504i
-0.21966219054509+0.56787488560383i
-0.21966219054509-0.56787488560383i
w=0.8时[x,n]=SOR(A,b,x0,0.8,10^(-5),300)
-7.99980126177984
0.333375649565234
3.66660553355675
1.99997665093436
238
w=1.2时
[x,n]=SOR(A,b,x0,1.2,10^(-5),300)
-7.99992226324999
0.333349197733832
3.66664330811682
1.99999119661784
111
w=1.3时
[x,n]=SOR(A,b,x0,1.3,10^(-5))
-7.99993630689315
0.333346074182271
3.66664773570359
1.99999290958378
89
w=1.6时
[x,n]=SOR(A,b,x0,1.6,10^(-5))
-7.99999410032529
0.333332406655426
3.66666433369261
1.99999911550516
26
比较上述方法可得,w=1.6时的SOR法最为快速。
(3)Cholesky分解求解下列线性方程组
A=[42-402400;
22-1-21320;
-4-1141-8-356;
0-216-1-4-33;
21-8-1224-10-3;
43-3-44111-4;
025-3-101142;
0063-3-4219];
b=[0;
-6;
20;
23;
9;
-22;
-15;
45];
Choleshy(A,b)
121.1481-140.112729.7515-60.152810.9120-26.79635.4259-2.0185
2.3设
(1)将矩阵A进行LU分解A=LU,其中U是上三角矩阵,L是主对角线上的元素都是1的下三角矩阵。
(2)利用上述分解分别求解方程组
并由此求出逆矩阵
。
(3)用LU分解求下列线性方程组的解
(方程组的精确解是
A=[2023;
181;
2-315];
[L,U,flag]=LU_Decom(A)
b1=[100]'
b2=[010]'
b3=[001]'
x1=A\b1,x2=A\b2,x3=A\b3
L=
1.000000
0.05001.00000
0.1000-0.40511.0000
U=
20.00002.00003.0000
07.90000.8500
0015.0443
OK
X=[x1,x2,x3]
X=
0.0517-0.0164-0.0093
-0.00550.1237-0.0072
-0.00800.02690.0665
X即为A的逆。
(2)matlab中输入命令:
[L,U,flag]=LU_Decom(A),formatlongg;
x=U\(L\b)得
L为
U为
1
-0.999999999999999
4.21884749357559e-015
0.999999999999999
1.99999999999999
1.68753899743024e-015
3
-1
2
2.4
(1)用追赶法求解方程组
(a)
(b)
(2)设计实验验证Hilbert矩阵
的病态性,其中
.
a1=4*ones(30,1)'
a2=ones(29,1)'
a3=a2;
b=ones(30,1)'
b(1,1)=2;
b(30,1)=2;
A=diag(a1)+diag(a3,1)+diag(a2,-1);
x=chasing(A,b)
zhuiganfa(A,b)
0.5359
-0.1436
0.0385
-0.0103
0.0028
-0.0007
0.0002
-0.0001
0.0000
-0.0000
0.5359
(b)
a1=5*ones(20,1)'
a2=-2*ones(19,1)'
a4=ones(18,1)'
a5=a4;
b=zeros(20,1);
b(1,1)=1;
A=diag(a1)+diag(a3,1)+diag(a2,-1)+diag(a4,-2)+diag(a5,2);
zhuiganfa(A,b)
0.249999999999829
0.124999999999574
0.0624999999991047
.0312********
0.0156249999963656
0.00781249999272582
0.00390624998544897
0.00195312497089661
0.000976562441792561
0.000488281133584789
0.000244140392169412
0.00012206984683874
6.10342249274393e-005
3.05157154798577e-005
1.5255063772205e-005
7.62194395065481e-006
3.79979610443202e-006
1.87754631042523e-006
8.94069671631063e-007
3.57627868652425e-007
matlab代码为:
m=input('
inputm:
='
%输入矩阵的阶数N=[m];
N=[m];
length(N)n=N(k);
%矩阵的阶
H=hilb(n);
%产生n阶Hilbert矩阵
disp(H)%输出n阶Hilbert矩阵
Hi=invhilb(n);
%产生完全准确的n阶逆Hilbert矩阵
b=ones(n,1);
%生成n阶全1向量
x_approx=H\b;
%利用左除H求近似解
x_exact=Hi*b;
%利用准确逆Hilbert矩阵求准确解
ndb=norm(H*x_approx-b);
nb=norm(b);
ndx=norm(x_approx-x_exact);
nx=norm(x_approx);
er_actual(k)=ndx/nx;
%实际相对误差
K=cond(H);
%计算Hilbert矩阵的条件数
er_approx(k)=K*eps;
%最大可能的近似相对误差
er_max(k)=K*ndb
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 数值 分析 课程设计 作业