数值分析课程设计Word文档下载推荐.docx
《数值分析课程设计Word文档下载推荐.docx》由会员分享，可在线阅读，更多相关《数值分析课程设计Word文档下载推荐.docx（25页珍藏版）》请在冰豆网上搜索。
首先用矩阵A的各阶主子式值判断矩阵能否进行LU分解，在能够分解的前提下，使用高斯消去法，将将消去时的两行间的倍数m附入L（m=A（k,p）/A（p,p）;
L（p,p）=1;
L（k,p）=m;
），再将消去后的A附入U，即得L,U。
function[L,U]gaussLU（A）
[nn]=size（A）;
RA=rank（A）;
ifRA~=n
disp（'
因为A的n阶行列式hl等于零，所以A不能进行LU分解。
'
）;
hl=det（A）;
return
ifRA==n
forp=1:
n
h（p）=det（A（1:
p,1:
p））;
hl=h（1:
n）;
fori=1:
ifh（1,i）==0
因为A的r阶主子式等于零，所以A不能进行LU分解。
）
ifh（1,i）~=0
因为A的各阶主子式都不等于零，所以A能进行LU分解。
L与U的值如下：
L=zeros（n,n）;
L（n,n）=1;
U=zeros（n,n）;
n-1
fork=p+1:
m=A（k,p）/A（p,p）;
L（p,p）=1;
A（k,p:
n）=A（k,p:
n）-m*A（p,p:
U=A（1:
n,1:
A=[2023;
181;
235];
[L,U]=gaussLU（A）
程序输出：
L=
1.000000
0.05001.00000
0.10000.35441.0000
U=
20.00002.00003.0000
07.90000.8500
004.3987
2．2用矩阵分解方法求上题中A的逆矩阵。
记
分别求解方程组
由于三个方程组系数矩阵相同，可以将分解后的矩阵重复使用。
对第一个方程组，由于A=LU，所以先求解下三角方程组
再求解上三角方程组
，则可得逆矩阵的第一列列向量；
类似可解第二、第三方程组，得逆矩阵的第二列列向量的第三列列向量。
由三个列向量拼装可得逆矩阵
。
我觉得题目中已经分析的很全面了，此题没有必要再多此一举了。
A=[20,2,3;
1,8,1;
2,-3,15];
[L,U]=gaussLU（A）;
b1=[1;
0;
0];
b2=[0;
1;
b3=[0;
1];
Y1=L\b1;
Y2=L\b2;
Y3=L\b3;
X1=U\Y1;
X2=U\Y2;
X3=U\Y3;
X=[X1,X2,X3]%X
输出结果：
X=
0.0517-0.0164-0.0093
-0.00550.1237-0.0072
-0.00800.02690.0665
实验三
3．1用泰勒级数的有限项逼近正弦函数
用计算机绘出上面四个函数的图形。
直接设定变量范围、函数，然后作图即可，注意图像可能会重叠，最好用不同线条，区别大一点的颜色绘制。
x=0:
pi/100:
pi;
x1=0:
pi/2;
y0=sin（x）;
y1=x1;
y2=x1-x1.^3/6;
y3=x1-x1.^3/6+x1.^5/120;
plot（x,y0,'
:
k'
x1,y1,'
-g'
x1,y2,'
-b'
x1,y3,'
-r'
输出图像：
3．2绘制飞机的降落曲线
一架飞机飞临北京国际机场上空时，其水平速度为540km/h，飞行高度为1000m。
飞机从距机场指挥塔的横向距离12000m处开始降落。
根据经验，一架水平飞行的飞机其降落曲线是一条三次曲线。
建立直角坐标系，设飞机着陆点为原点O，降落的飞机为动点
，则
表示飞机距指挥塔的距离，
表示飞机的飞行高度，降落曲线为
该函数满足条件：
（1）试利用
满足的条件确定三次多项式中的四个系数；
（2）用所求出的三次多项式函数绘制出飞机降落曲线。
通过y（x）满足的条件带入降落曲线方程，即可得到一个四元一次方程组。
将方程组视为矩阵。
由于前面的题目用过了高斯消去法，本体继续沿用此方法，通过高斯消去法，求解方程组，得到a0,a1,a2,a3的值，带回y中，运用plot函数，输出x，y的关系曲线。
functiona=jiangl（x1,x2,x3,x4,b）
A=[1x1x1^2x1^3;
1x2x2^2x2^3;
012*x33*x3*x3;
012*x43*x4*x4];
B=[Ab];
a=zeros（4,1）;
C=zeros（1,5）;
forp=1:
3
[Y,j]=max（abs（B（p:
4,p）））;
C=B（p,:
B（p,:
）=B（j+p-1,:
B（j+p-1,:
）=C;
4
m=B（k,p）/B（p,p）;
B（k,p:
5）=B（k,p:
5）-m*B（p,p:
5）;
b=B（1:
4,5）;
A=B（1:
4,1:
4）;
a（4）=b（4）/A（4,4）;
forq=3:
-1:
1
a（q）=（b（q）-sum（A（q,q+1:
4）*a（q+1:
4）））/A（q,q）;
在指令窗口输入：
x1=0;
x2=12000;
x3=0;
x4=12000;
b=[0;
1000;
a=jiangl（x1,x2,x3,x4,b）
a=
1.0e-004*
0
0.2083
-0.0000
x=0:
50:
12000;
y=a
（1）+x*a
（2）+x.^2*a（3）+x.^3*a（4）;
plot（x,y）
回车，输出图像：
4．1曾任英特尔公司董事长的摩尔先生早在1965年时，就观察到一件很有趣的现象：
集成电路上可容纳的零件数量，每隔一年半左右就会增长一倍，性能也提升一倍。
因而发表论文，提出了大大有名的摩尔定律（Moore’sLaw），并预测未来这种增长仍会延续下去。
下面数据中，第二行数据为晶片上晶体数目在不同年代与1959年时数目比较的倍数。
这些数据是推出摩尔定律的依据：
年代
1959
1962
1963
1964
1965
增加倍数
5
6
试从表中数据出发，推导线性拟合的函数表达式。
通过表中数据绘图，可发现趋向于一次线性方程，正好本题要求推到线性拟合函数，直接使用一次的最小二乘法。
functionzxec（x,y）
y=y'
;
fori=1:
length（x）
forj=1:
2
A（i,j）=x（i）^（j-1）;
[L,U]=lu（A'
*A）;
alpha=U\（L\（A'
*y））
x=[1959,1962,1963,1964,1965];
y=[1,3,4,5,6];
zxec（x,y）
alpha=
1.0e+003*
-1.6255
0.0008
其中alpha
（1）=-1625.5是函数交y轴的焦点值，0.8是一阶函数斜率
5．1用几种不同的方法求积分
的值。
（1）牛顿-莱布尼茨公式；
（2）梯形公式；
（3）辛卜生公式；
（4）复合梯形公式。
（1）牛顿-莱布尼茨
由于
=4arctgx|0到1
输入内容：
4*atan
（1）
ans=
3.1416
（2）（3）（4）
定义
functions=f（x）
s=4/（1+x^2）;
输入：
a=0;
b=1*
（2）梯形公式
接着*式输入：
（1/2）*（b-a）*（f（b）+f（a））
3
（3）辛卜生公式
（1/6）*（b-a）*（f（a）+4*f（（a+b）/2）+f（b））
3.1333
（4）复合梯形公式
T=0;
forn=0:
0.1:
0.9
m=0.1/2*（f（n）+f（n+0.1））;
T=T+m;
T
T=
3.1399
6．1用欧拉公式和四阶龙格-库塔法分别求解下列初值问题；
此题只需建立欧拉和龙格库塔函数，带入求解即可。
由于题目没有给出步长，我们取步长h=0.1
（1）建立名为f.m的文件：
functionr=f（t,y）
r=（0.9*y）/（1+2*t）;
Euler方法
functioneuler（a,b,h,alpha）
t=a:
h:
b;
n=（b-a）/h;
y
（1）=alpha;
fori=2:
n+1
y（i）=y（i-1）+h*f（t（i-1）,y（i-1））;
sprintf（'
t=%3.1f,y=%9.7f'
t（i）,y（i））
在命令窗口输入：
euler（0,1,0.1,1）
回车得到：
t=0.1,y=1.0900000
t=0.2,y=1.1717500
t=0.3,y=1.2470768
t=0.4,y=1.3172249
t=0.5,y=1.3830861
t=0.6,y=1.4453250
t=0.7,y=1.5044519
t=0.8,y=1.5608688
t=0.9,y=1.6148989
t=1.0,y=1.6668064
四阶Runge-Kutta方法
functionrk4（a,b,h,alpha）
k1=f（t（i-1）,y（i-1））;
k2=f（t（i-1）+h/2,y（i-1）+h/2*k1）;
k3=f（t（i-1）+h/2,y（i-1）+h/2*k2）;
k4=f（t（i-1）+h,y（i-1）+h*k3）;
y（i）=y（i-1）+h/6*（k1+2*k2+2*k3+k4）;
rk4（0,1,0.1,1）
t=0.1,y=1.0855051
t=0.2,y=1.1634777
t=0.3,y=1.2355334
t=0.4,y=1.3027862
t=0.5,y=1.3660420
t=0.6,y=1.4259056
t=0.7,y=1.4828446
t=0.8,y=1.5372290
t=0.9,y=1.5893579
t=1.0,y=1.6394763
同理：
（2）建立名为f.m的文件：
r=-（t*y）/（1+t^2）;
在命令窗口分别输入：
euler（0,1,0.1,2）
回车
rk4（0,1,0.1,2）
结果如下表
t=0.1,y=2.0000000
t=0.2,y=1.9801980
t=0.3,y=1.9421173
t=0.4,y=1.8886645
t=0.5,y=1.8235382
t=0.6,y=1.7505966
t=0.7,y=1.6733644
t=0.8,y=1.5947500
t=0.9,y=1.5169573
t=1.0,y=1.4415285
t=0.1,y=1.9900743
t=0.2,y=1.9611612
t=0.3,y=1.9156523
t=0.4,y=1.8569530
t=0.5,y=1.7888539
t=0.6,y=1.7149854
t=0.7,y=1.6384634
t=0.8,y=1.5617372
t=0.9,y=1.4865879
t=1.0,y=1.4142132
7.1用二分法求下列方程在指定区间[a，b]上的实根近似值：
（要求误差不超过0.01）
（1）x-sinx-1=0，[a,b]=[1,3]；
（2）xsinx=1，[a,b]=[0,1.5].
直接将方程带入二分法函数中求解即可。
functionx=eff（fun,a,b,eps）
ifnargin<
eps=1e-5;
fa=feval（fun,a）;
fb=feval（fun,b）;
iffa*fb>
[a,b]不是有根区间，请重新调整'
return;
whileabs（b-a）/2>
eps
x=（a+b）/2;
fx=feval（fun,x）;
iffx*fa<
b=x;
fb=fx;
else
a=x;
fa=fx;
x=（a+b）/2;
（1）fun=inline（'
x-sin（x）-1'
x=eff（fun,1,3,0.01）
输出结果
x=
1.9297
（2）fun=inline（'
x*sin（x）-1'
x=eff（fun,0,1.5,0.01）
1.1191
8.1分别用直接法、雅可比迭代法、赛德尔迭代法求解线性方程组AX=b。
（1）
，
（2）
首先通过diag函数给矩阵A附上初值，然后分别调用雅可比和赛德尔函数即可
程序：
M=ones（9,9）;
N=4*ones（10,10）;
A=diag（diag（N））+diag（diag（M）,1）+diag（diag（M）,-1）;
b=[2;
2];
1.用直接法输入：
formatlong
A\b
0.479271991911021
0.082912032355915
0.189********5319
0.160768452982811
0.167846309403438
2.用Jacobi法
functionJacobi（A,b,X0,P,error,maxl）
X=zeros（n,1）;
fork=1:
maxl
X（j）=（b（j）-A（j,[1:
j-1,j+1:
n]）*X0（[1:
n]））/A（j,j）;
errX=norm（X-X0,P）;
X0=X;
if（errX<
error）
迭代次数k，近似解X分别是：
k
X
if（errX>
=error）
请注意：
Jacobi迭代次数已经超过最大迭代次数maxl.'
X0=[0;
Jacobi（A,b,X0,inf,0.001,100）;
迭代次数k，精确解X1和近似解X分别是：
k=
9
0.47933959960938
.0831********
0.189********852
0.16108322143555
0.16813659667969
用赛德尔迭代法：
functionGS（A,b,X0,P,error,maxl）
XX=0;
ifi<
j
XX=XX+A（j,i）*X（i）;
ifi>
XX=XX+A（j,i）*X0（i）;
X（j）=（b（j）-XX）/A（j,j）;
Gauss-Seidel迭代次数已经超过最大迭代次数maxl.'
GS（A,b,X0,inf,0.001,100）;
6
0.47925853729248
0.08289575576782
0.189********913
0.16064277291298
0.16802297346294
0.16770140314475
0.16085489129182
0.189********114
0.08292683955369
0.47926829011158
U=ones（18,18）;
V=-2*ones（19,19）;
W=5*ones（20,20）;
A=diag（diag（W））+diag（diag（V）,1）+diag（diag（U）,2）+diag（diag（V）,-1）+diag（diag（U）,-2）;
b=[1;
zeros（19,1）];
X0=zeros（20,1）;
同上获得如下解：
直接法
雅可比
赛德尔
7
8
9
10
11
12
13
14
15
16
17
18
19
20
0.242121373548085
.0943********
.021*********
.0313********
-0.006942557862112
0.004886927715767
0.003586117214438
0.000214823135181
-0.000784526508185
-0.000357861979163
0.000050437638520
0.000106309021306
0.000029232399870
-0.000014343050585
-0.000012667696430
-0.000001464361499
0.000002491258113
0.000001311981447
-0.000000093354238
-0.000000299737984
1.0e+004*-0.28210827857172
1.0e+004*0.48204640630127
1.0e+004*-0.69195375595679
1.0e+004*0.88126694207517
1.0e+004*-1.05291766695873
1.0e+004*1.20130020954019
1.0e+004*-1.32375170735604
1.0e+004*1.41753761486760
1.0e+004*-1.48072923013544
1.0e+004*1.51202411408633
1.0e+004*-1.51082661747071
1.0e+004*1.47723879582980
1.0e+004*-1.41205236153589
1.0e+004*1.31674011324409
1.0e+004*-1.19335996640489
1.0e+004*1.04472588730159
1.0e+004*-0.87353157471502
1.0e+004*0.68530915502268
1.0e+004*-0.47711982856112
1.0e+004*0.27913406015473
0.24229171609600
0.09477250842624
.0314********
-0.0071474575114
升级会员 