数值分析课程课程设计.docx

数值分析课程课程设计.docx

《数值分析课程课程设计.docx》由会员分享，可在线阅读，更多相关《数值分析课程课程设计.docx（24页珍藏版）》请在冰豆网上搜索。

数值分析课程课程设计

计目设题

共15题如下

成绩

课程设计

我再也回不到大二了，

大学是那么短暂

设计题目数值分析

学生姓名李飞吾

字亏xxxxxxxx

专业班级信息计xxxxx班

指导教师

课程设计主要内容

指导老师评语

设计目的：

通过不同题目的理解，进行算法分析。

通过MATLAB软件进行编程对题目进行解决。

个别题目设计验证，加深对数值分析的理解。

函数的图像绘制的运用

设计题目：

1题1.1利用逆向递推的方法求解问题，通过条件终止地推

2题1.2从某个初始值开始，利用递推公式进行积分估值

@题1.3绘制Koch分形曲线，节点之间的关系与坐标变换

（4题2.1用高斯消元法的消元过程作矩阵分解，LU分解

⑤题2.2矩阵分解方法求上题中A的逆矩阵，针对不同的b,而重复利用已知的LU

®题2.3验证希尔伯特矩阵的病态性，矩阵基本运算

7题3.1用泰勒级数的有限项逼近正弦函数，由图像观察逼近效果

8题3.2绘制飞机的降落曲线，线性方程组求解，与绘图

⑤题4.1线性拟合的函数表达式的推导，使用了两种代码方法

面题5.1用几种不同的方法求积分，观察数值积分的逼近效果

®题5.5求空间曲线L弧长。

求导后使用符号函数积分计算

O题6.1用欧拉公式和四阶龙格-库塔法分别求解下列初值问题，代码搜索内容。

（0题6.4常微分方程的解，dsolve（）函数使用

@）题8.2差分法解常微分方程边值问题，ode函数无能为力，Matlab中提供bvp

食军算器。

solinit=bvpinit（x,yinit,params）sol=bvpsolver（odefun,bcfun,solinit,options）

（g题8.3求解圆的半径，圆心。

线性方程组解参数

设计总结：

（1）算法是题目的解题核心，好的算法可以使计算更加精确（题5.1）

（2）图形绘制在今后的课程设计，或者是论文中可以用到。

（3）无法解决的问题，需要请教室友，或者上网查阅。

（4）MATLAB是一个很强大的软件，提供了很多内置的数学函数，直接进行解题。

查阅资料时了解到一些MATLAB论坛。

通过帖子阅读，了解到了MATLAB在科

学计算方面，模拟，图形，视频等起到的作用。

增加了对“计算科学“的理解。

建议：

从学生的工作态度、工作量、设计（论文的）创造性、学术性、使用性及书面表达能力等方面给出评价。

签名：

20年月日

数值分析谋程设计

1.1水手、猴子和椰子问题：

五个水手带了一只猴子来到南太平洋的

一个荒岛上，发现那里有一大堆椰子。

由于旅途的颠簸，大家都很疲惫，很快就入睡了。

第一个水手醒来后，把椰子平分成五堆，将多余的一只给了猴子，他私藏了一堆后便又去睡了。

第二、第三、第四、第五个水手也陆续起来，和第一个水手一样，把椰子分成五堆，恰多一只猴子，私藏一堆，再去入睡，天亮以后，大家把余下的椰子重新等分成五堆，每人分一堆，正好余一只再给猴子，试问原先共有几只椰子？

（15621）

试分析椰子数目的变化规律，利用逆向递推的方法求解这一问题解：

算法分析：

解该问题主要使用递推算法，关于椰子数目的变化规律可以设起初的椰子数为po,第一至五次猴子在夜里藏椰子后，椰子的数目分别为po,p1,p2,P3,P4再设最后每个人分得x个椰子，由题：

r17（pk1）（k=0,1,2,3,4）x-（p51）

55

所以P55x1,p

R5妇1,（k=0,1,2,3,4）

4

输出结果

•

1023

15621

MATLAB码：

n=input（'n='）;n=15621forx=1:

np=5*x+1;fork=1:

5p=5*p/4+1;

end

结果分析：

此题的解题思想是在迭代法中，判断p为整数时，输出与p

ifp==fix（p）,breakendenddisp（[x,p]）

dn

1x

1.2设，Indx

05x

（1）从10尽可能精确的近似值出发，利用递推公式：

1In5In1—（n1,2,L20）n

计算机从Ii到I20的近似值；

（2）从I30较粗糙的估计值出发，用递推公式：

Ini-In—（n30,29,L,3,2）

55n

计算从Il到I20的近似值；

解：

首先第一步，估计Io和I30的值：

symsxn;

int（xA0/（5+x）,0,1）

ans=log

（2）+log（3）-log（5）

eval（ans）

ans=

0.1823

则取I。

为0.18

symsxn;

int（xA30/（5+x）,0,1）ans=

931322574615478515625*log

（2）+931322574615478515625*log（3）-931

322574615478515625*log（5）-79095966183067699902965545527073/46

5817912560

eval（ans）

ans=

0

MATLAB中小数点后保留四位，由上面计算知道积分的值不为了零。

所以I30的取值为0.00001-0.0001

MATLAB码：

i=input（'i='）;

i=0.18;

>>ifi>=0.1&&i<=0.2

forn=1:

1:

20

i=-5*i+1/n

end

elseifi>0&&i<=0.0001

forn=30:

-1:

2

i=（-1/5）*i+1/（5*n）

end

end

i=

1.1336e+005

i=

-5.6679e+005

i=

2.8339e+006

i=

-1.4170e+007

i=

7.0848e+007

i=

-3.5424e+008

i=

1.7712e+009

i=

-8.8560e+009

i=

4.4280e+010

i=

-2.2140e+011

同理输入积分初始值i=0时可以得i=0.0884

结果分析：

第二种方法所得的结果相对来说比较精确一些，也比较可靠因为第一种方法每一迭代都将最初的误差放大了五倍，使得最终的误差

越来越大；而第一种方法经过每一次迭代都将误差缩小为初始误差的五分之一，使得最终的误差越来越小，因此相对来说比较可靠，性能较好。

1.3绘制Koch分形曲线

问题描述：

从一条直线段开始，将线段中间的三分之一部分用一个等边

三角形的另两条边代替，形成具有5个结点C0SJ的国写（图1-4）;在新的图形中，又将图中每一直线段中间的询之丁部分都用一个等边三角形的另两条边代替，再次形成新的图形（图sit5）,cos时，图形中共有17

个结点。

这种迭代继续进行下去可以形成Koch分形曲线。

问题分析：

考虑由直线段（2个点）产生第一个图形（5个点）的过程，设p1和R分别为原始直线段的两个端点。

现在需要在直线段的中间依次插入三个点P2，R，P4产生第一次迭代的图形（图1-4）。

显然，只位于P点右端直线段的三分之一处，已点绕p2旋转60度（逆时针方向）而得到的，故可以处理为向量巳已经正交变换而得到向量P2P3,形成算法如下：

（1）

P2

P（任P）/3.

（2）

R

P2（P5R）/3.

（3）

P3

%（EE）AT.

在算法的第三步中，A为正交矩阵。

这一算法将根据初始数据（P和R点的坐标），产生图1-4中5个结点的

坐标。

这5个结点的坐标数组，组成一个5X2矩阵。

这一矩阵的第一行为为R的坐标，第二行为R的坐标，第二行为P2的坐标……第五行为P5的坐标。

矩阵的第一列元素分别为5个结点的x坐标，第二列元素分别为5个结点的y坐标。

问题思考与实验：

（1）考虑在Koch分形曲线的形成过程中结点数目的变化规律。

设第k次迭代产生结点数为nk,第k1迭代产生结点数为nk1,试写出nk和nk1之间的递推关系式；

（2）参考问题分析中的算法，考虑图1-4到图1-5的过程，即由第一次迭代的5个结点的结点坐标数组，产生第二次迭代的17个结点的结点坐标数组的算法；

（3）考虑由第k次迭代的nk个结点的结点坐标数组，产生第k1次迭代的nk1个结点的结点坐标数组的算法；

（4）设计算法用计算机绘制出如下的

解：

（1）nk14nk3

Koch分形曲线（图1-6）

（2）（3辩法及（4）代码分析：

p=[00;100];%P为初始两个点的坐标，第一列为x坐标，第二列为y坐标

n=2;%n为结点数

A=[cos（pi/3）-sin（pi/3）;sin（pi/3）cos（pi/3）];%旋转矩阵

fork=1:

4

d=diff（p）/3;%diff%

m=4*n-3;%

q=p（1:

n-1,:

）;%

p（5:

4:

m,:

）=p（2:

n,:

）;%p（2:

4:

m,:

）=q+d;%

p（3:

4:

m,:

）=q+d+d*A';%p（4:

4:

m,:

）=q+2*d;%n=m;%

end

plot（p（:

1）,p（:

2））%

axis（[01003.5]）

计算相邻两个点的坐标之差，得到相邻两点确定的向量

则d就计算出每个向量长度的三分之一，与题中将线段三等分对应

迭代公式

以原点为起点，前n-1个点的坐标为终点形成向量

迭代后处于4k+1位置上的点的坐标为迭代前的相应坐标

用向量方法计算迭代后处于

用向量方法计算迭代后处于

用向量方法计算迭代后处于

4k+2位置上的点的坐标

4k+3位置上的点的坐标

4k位置上的点的坐标

迭代后新的结点数目

2.1用高斯消元法的消元过程作矩阵分解。

设

2023

A181

2315

消元过程可将矩阵A化为上三角矩阵U,试求出消元过程所用的乘数

m21、m31、伉1并以如下格式构造下三角矩阵L和上三角矩阵U

12023

Lm211,Ua22）a23）

4

（2）

m31m321a33

验证：

矩阵A可以分解为L和U的乘积，即A=LU。

矩阵LU分解MATLAB代码：

functionhl=zhjLU（A）

[n,n]=size（A）;RA=rank（A）;

ifRA~=n

disp（'因为A的n阶行列式hl等于零，所以A不能进行LU分解.A的

秩RA如下：

’）；

RA,hl=det（A）;

return

end

ifRA==n

forp=1：

n

h（p）=det（A（1：

p,1：

p））;

end

hl=h（1：

n）;

fori=1：

n

ifh（1,i）==0

disp（'因为A的各阶主子式等不等于零，所以A能进彳fLU

分解.A的秩RA和各阶顺序主子式值hl庆次如下：

’）；

RA,hl

return

end

end

ifh（1,i）~=0

disp（'因为A的各阶主子式都不等于零，所以A能进行LU分解.A的秩RA和各阶顺序主子式值hl如下：

’）；

forj=1：

n

U（1,j）=A（1,j）;

end

fork=2:

n

fori=2:

n

forj=2:

n

L（1,1）=1;L（i,i）=1;ifi>j

L（1,1）=1;L（2,1）=A（2,1）/U（1,1）;L（i,1）=A（i,1）/U（1,1）;

L（i,k）=（A（i,k）-L（i,1:

k-1）*U（1:

k-1,k））/U（k,k）;

else

U（k,j）=A（k,j）-L（k,1:

k-1）*U（1:

k-1,j）;

end

end

end

end

RA,hl,U,L

end

end

仅上上代码保存为M文件，并在命令窗口输入

A=[2023;181;2-315];

b=[000]';

h=zhjLU（A）

程序运行结果:

L=

U=

1.0000

0020.00002.00003.0000

0.0500

1.0000007.90000.8500

0.1000

-0.40511.00000015.0443

实验验证：

可以直接使用MATLAB内置LU分解

A=[2023;181;2-315];

[LU]=lu（A）输出结果与上程序输出结果一致。

2.2用矩阵分解方法求上题中A的逆矩阵。

记

1

0

0

b0,b2

1,b3

0

0

0

1分别求解方程组AXb1，AXb2，AX

由于三个方程组系数矩阵相同，可以将分解后的矩阵重复使用。

对第一

个方程组，由于A=LU，所以先求解下三角方程组LY知再求解上三角方程组UXY,则可得逆矩阵的第一列列向量；类似可解第二、第三方程组，得逆矩阵的第二列列向量的第三列列向量。

由三个列向量拼装可得逆矩阵A1。

解:

MATLAB代码如下:

b1=[1;0;0];b2=[0;1;0];b3=[0;1;1];

A=[20,2,3;1,8,1;2,-3,15];

L=[1,0,0;0.05,1,0;0.1,-0.4051,1];

U=[2023;07.90.85;0015.0443];

Y1=L\b1

X1=U\Y1

Y2=L\b2X2=U\Y2Y3=L\b3

X3=U\Y3

Y1=

1.0000

-0.0500

-0.1203

X1=

0.0517

-0.0055

-0.0080

Y2=

0

实验验证：

[X1X2X3]ans=

0.0517-0.0164-0.0257

-0.00550.12370.1165

-0.00800.02690.0934

而：

inv（A）=[X1X2X3]得证

1.0000

0.4051

X2=

-0.0164

0.1237

0.0269

Y3=

0

1.0000

1.4051

X3=-0.02570.11650.0934

2.3验证希尔伯特矩阵的病态性：

对于三阶矩阵

11/21/3

H1/21/31/4

1/31/41/5

取右端向量b【11/613/1247/6°]t,验证：

（1）向量X[X1X2X3]T[111]T是方程组HXb的准确解；

（2）取右端向量b的三位有效数字得b[1.83L080.783]T,求方程组的准确解X,并与X的数据I1,1,1】'作比较。

说明矩阵的病态性。

解：

⑴

H=[11/21/3;1/21/31/4;1/31/41/5];

X=[1；1；1]；

b=H*X

b=

1.8333

1.0833

0.7833

与题中相同

（2）先求出解X,与数据[〔，"亍作比较

H=[11/21/3;1/21/31/4;1/31/41/5];

b=[1.83;1.08;0.783];

X=H\b

X=

1.0800

0.5400

1.4400

与[1,1,1]相差较大，矩阵为病态矩阵

3.

1用泰勒级数的有限项逼近正弦函数

sinx,x[0,]

x,x[0,/2]

3——

xx/6,x[0,/2]xx3/6x5/120,x

用计算机绘出上面四个函数的图形。

解:

MATLAB代码如下

（1）symsx;taylor（sin（x））x=0:

0.01*pi:

piplot（x,sin（x））

⑵

symsx;

taylor（x）

x=0:

0.01*pi:

pi/2

plot（x）

⑶

symsx;

taylor（x-xA3/6）

fplot（'x-xA3/6',[0pi/2]）（4）

symsx;

taylor（x-xA3/6+xA5/120）

fplot（'x-xA3/6+xA5/120',[0pi/2]）

结果图形右：

x=0:

0.1:

pi;

y=sin（x）;

plot（x,y,'-sk'）;

holdon

x=0:

0.1:

pi/2;

y=x;

plot（x,y,'-b*'）

holdon

fplot（'x-x.A3/6',[0,pi/2,0,2],2e-3,'-gx'）

holdonfplot（'x-x.A3/6+x.A5/120',[0,pi/2,0,2],2e-3,'-r.'）holdoff

legend（'sin（x）','x','x-xA3/6','x-A3/6+xA5/120',2）

xlabel（'x'）

ylabel（'y'）

title（'Taylorapproximation'）

结果分析：

图中红色点线为正弦曲线，蓝色的星线为一阶泰勒逼近，

绿色叉线为二阶

豪勒逼近，黑色正方形线为三阶泰勒逼近。

可见三阶泰勒逼近效果最好，泰勒级改越高，逼近效果尬婷。

3.2绘制飞机的降落曲线

一架飞机飞临北京国际机场上空时，其水平速度为540km/h,飞行高度为1000m。

飞机从距机场指挥塔的横向距离12000m处开始降落。

根据经验，一架水平飞行的飞机其降落曲线是一条三次曲线。

建立直角坐标系，设飞机着陆点为原点O，降落的飞机为动点P（x,y），则x表示飞机距指挥塔的距离，y表示飞机的飞行高度，降落曲线为

23

y（x）a0a1xa2xa3x

该函数满足条件：

y（0）0,y（12000）1000

y'（0）0,y'（12000）0

（1）试利用y（x）满足的条件确定三次多项式中的四个系数；

（2）用所求出的三次多项式函数绘制出飞机降落曲线。

functionS=f（x1,y1,x2,y2,x3,y3,x4,y4）formatlong

a1=[1,x1,x1A2,x1A3];

a2=[1,x2,x2A2,x2A3];

a3=[0,1,2*x3,3*x3A2];

a4=[0,1,2*x4,3*x4A2];

a=[a1;a2;a3;a4];

b=[y1;y2;y3;y4];

s=a\b;

x=-12000:

250:

0;

y=s（3）*x.A2-s（4）*x.A3

plot（x,y,'-d'）

xlabel（'x'）

ylabel（'y'）

以上为M文件内容，在命令窗口键入

f（0,0,12000,1000,0,0,12000,0）

运行结桌：

ans为多项土索数

ans=

1.0e-004*

0

0

0.20833333333333

-0.00001157407407

4.1曾任英特尔公司董事长的摩尔先生早在1965年时，就观察到一件很有趣的现象：

集成电路上可容纳的零件数量，每隔一年半左右就会增长一倍，性能也提升一倍。

因而发表论文，提出了大大有名的

年代

1959

1962

1963

1964

1965

增加倍数

1

3

4

5

6

摩尔定律（Moore'sLaW,并预测未来这种增长仍会延续下去。

下面数据中，第二行数据为晶片上晶体数目在不同年代与1959年时数目

比较的倍数。

这些数据是推出摩尔定律的依据：

试从表中数据出发，推导线性拟合的函数表达式。

解：

解法一

MATLAB：

码：

x=[1959,1962,1963,1964,1965];

y=[1,3,4,5,6];

p1=polyfit（x,y,1）

y1=polyval（p1,x）

plot（x,y1,'-',x,y,'r*'）

xlabel（'x'）,ylabel（'y'）;

运行结果：

p1=

1.0e+003*

0.0008-1.6255

y1=

0.81133.30194.13214.96235.7925

线性拟合的函数表达式：

Y=0.8302x-1.6255e+003

解法二：

x=[19591962196319641965];

运行结果：

xs=

1.0e+003*

-1.625528301891238

0.000830188679248

y=[1;3;4;5;6];

fori=1:

length（x）

forj=1:

2

A（i,j）=x（i）”1）;

end

end

[L,U]=lu（A'*A）;

xs=U\（L\（A'*y））

从而年代y与增加倍数x之间的关系为：

y=-1625.528301891238+0.830188679248x

运行结果：

11=

pi

12=3

13=

3.1333

14=

3.1416

（1）牛顿-莱布尼茨公式；

（2）梯形公式；（3）辛卜生公式；（4）复合梯形公式。

解：

symsx

i1=int（4/（1+xA2）,x,0,1）

a=0;b=1;

h=b-a;

i2=（4/（1+aA2）+4/（1+bA2））/2

i3=h/6*（4/（1+aA2）+4*4/（1+（a+b/2）A2）+4/（1+bA2））

M=100;

h=（b-a）/M;

i4=0;

fork=1:

（M-1）

x=a+h*k;

i4=i4+4/（1+x"2）;

end

i4=h*（4/（1+aA2）+4/（1+bA2））/2+h*i4

牛顿莱布尼兹公式得到精确结果a采用挣形公式得到的结果比采用Simpson公式的精蒲I度鱼很，米用复祝榜形公式在步长取得越束裁小的状态下可以提高精度。

5.5求空间曲线L:

ysint

z2costsint

弧长公式为

L:

.x'2（t）y'2（t）z'2（t）dt

xcost

6.1用欧拉公式和四阶龙格-库塔法分别求解下列初值问题;

0.9y

（i）y',y（0）i;x[0,1]

12x

2；x[0,1]

⑵y'闩顶。

）

1x

解：

（1）<1>欧拉公式：

function[t,x]=Euler（fun,t0,tt,x0,N）h=（tt-t0）/N;

t=t0+[0:

N]'*h;

x（1,:

）=x0';

fork=1:

N

f=feval（fun,t（k）,x（k,:

））;

f=f;

x（k+1,:

）=x（k,:

）+h*f;

end

以'Euler.m'保存

functionf=Euler_fun（t,x）f=0.9*x./（1+2*t）「以'Euler_fun.m'保存functionmain_Euler[t,x]=Euler（'Euier_fun',0,1,1,20）;

fh=dsolve（'Dx=0.9*x/（1+2*t）','x（0）=1'）;fork=1:

8

ft（k）=t（k）;

fx（k）=subs（fh,ft（k））;

end

[t,x],4

以'main_Euler.m'保存输入：

main_Euler

<2>四阶龙杯-库塔法

functionR=rk4（f,a,b,ya,N）

%y'=f（x,y）

%a,b为左右端点

%Nf