数学分析课程设计1.docx
《数学分析课程设计1.docx》由会员分享，可在线阅读，更多相关《数学分析课程设计1.docx（69页珍藏版）》请在冰豆网上搜索。
实验一1.1 水手、猴子和椰子问题：
五个水手带了一只猴子来到南太平洋的一个荒岛上，发现那里有一大堆椰子。
由于旅途的颠簸，大家都很疲惫，很快就入睡了。
第一个水手醒来后，把椰子平分成五堆，将多余的一只给了猴子，他私藏了一堆后便又去睡了。
第二、第三、第四、第五个水手也陆续起来，和第一个水手一样，把椰子分成五堆，恰多一只猴子，私藏一堆，再去入睡，天亮以后，大家把余下的椰子重新等分成五堆，每人分一堆，正好余一只再给猴子，试问原先共有几只椰子？
试分析椰子数目的变化规律，利用逆向递推的方法求解这一问题。
解：
问题分析：
假设第n个人分椰子前，有An个椰子，即第n-1醒来面对的椰子数为An-1，第n个水手醒来后，把椰子平分成五堆，将多余的一只给了猴子，他私藏了一堆后便又去睡了则易得：
（1）
设最后一堆为a个，则第五个人分完后余下A（6）=5*a+1;由
（1）式可得出A
（1）,A
（2）...A（5）;
由于椰子数目为整数，有a=1开始，建立循环，可以由matlab求得最小的a=1023 ，五次数量变化如下：
Matlab窗口输入如下：
>>a=1;
A=zeros（1,6）;
whilea>0,
A（6）=5*a+1;
if A（6）~=fix（A（6））,
a=a+1;
continue;
end
i=6;
while i>=2,
A（i-1）=A（i）*5/4+1;
if A（i-1）~=fix（A（i-1））,
break;
end
i=i-1;
end
if （A
（1）~=0）&&（A
（1）==fix（A
（1）））,
break;
end
a=a+1;
end
>>a
a=
1023
>>A
（1）
ans
15621
>>plot（A,’k’）
1.2 当时，选择稳定的算法计算积分
解：
用蒙特卡罗方法求解该问题，在（0，1）上随机选取一定数量的点，以落在积分区域内的频率近似所求积分，matlab程序如下：
mtj.m
functiona=mtj（n）
f=@（x）exp（-n*x）/（exp（-x）+10）;
N=10000;
r=0;
fori=1:
N,
A=rand（1,2）;
if A
（2）<=f（A
（1））
r=r+1;
end
end
a=r/N;
1.3 绘制静态和动态的Koch分形曲线
问题描述：
从一条直线段开始，将线段中间的三分之一部分用一个等边三角形的另两条边代替，形成具有5个结点的新的图形；在新的图形中，又将图中每一直线段中间的三分之一部分都用一个等边三角形的另两条边代替，再次形成新的图形，这时，图形中共有17个结点。
这种迭代继续进行下去可以形成Koch分形曲线。
在迭代过程中，图形中的结点将越来越多，而曲线最终显示细节的多少取决于所进行的迭代次数和显示系统的分辨率。
Koch分形曲线的绘制与算法设计和计算机实现相关。
图1.1 Koch曲线的形成过程
解：
A.绘制静态Koch分形曲线，代码如下：
A=[cos（pi/6）-sin（pi/6）;sin（pi/6） cos（pi/6）];
p=[010;00;];%第一行是横坐标第二行是纵坐标
n=2;
fork=1:
6,
m=1;
fori=1:
n-1,
d=p（:
i+1）-p（:
i）;
p1（:
m+1）=p（:
i）+d/3;
p1（:
m+3）=p（:
i）+d*2/3;
p1（:
m+2）=p（:
i）+A*d/sqrt（3）;
p1（:
m:
4:
m+4）=p（:
i:
i+1）;
m=m+5;
end
n=length（p1）;
p=p1;
clearp1;
end
x=p（1,:
）;
y=p（2,:
）;
plot（x,y,'k'）;
axis（[010010]）;
matalb得到图形如下：
B.动态的Koch分形曲线，在静态曲线的基础上得到动态分形曲线的代码和截图如下：
A=[cos（pi/6）-sin（pi/6）;sin（pi/6）cos（pi/6）];
p=[010;00;];
n=2;
fork=1:
6,
m=1;
fori=1:
n-1,
d=p（:
i+1）-p（:
i）;
p1（:
m+1）=p（:
i）+d/3;
p1（:
m+3）=p（:
i）+d*2/3;
p1（:
m+2）=p（:
i）+A*d/sqrt（3）;
p1（:
m:
4:
m+4）=p（:
i:
i+1）;
m=m+5;
end
n=length（p1）;
p=p1;
clearp1;
x=p（1,:
）;
y=p（2,:
）;
plot（x,y,'k'）;
axis（[010010]）;
pause
（2）;
end
实验二2.1 小行星轨道问题：
一天文学家要确定一颗小行星绕太阳运行的轨道，他在轨道平面内建立以太阳为原点的直角坐标系，在五个不同的对小行星作了五次观察，测得轨道上五个点的坐标数据（单位：
万公里）如下表所示：
P1
P2
P3
P4
P5
X坐标
53605
58460
62859
66662
68894
Y坐标
6026
11179
16954
23492
68894
由开普勒第一定律知，小行星轨道为一椭圆，椭圆的一般方程可表示为：
现需要建立椭圆的方程以供研究。
（1） 分别将五个点的数据代入椭圆一般方程中，写出五个待定系数满足的等式，整理后写出线性方程组
以及方程组的系数矩阵和右端项b；
（2） 用MARLAB求低阶方程的指令A＼b求出待定系数;
（3） 分别用直接法、Jacobi迭代法、Gauss-Seidel迭代法求出待定系数.
解：
（1）由五个点的坐标（x 1,y1）...（x5,y5）代入一般椭圆方程整理可得
代入题目中的数据得
A=[2873496025.0, 646047460.0, 36312676.0,107210.0, 12052.0]
[3417571600.0,1307048680.0, 124970041.0,116920.0, 22358.0]
[3951253881.0,2131422972.0, 287438116.0,125718.0, 33908.0]
[4443822244.0,3132047408.0, 551874064.0,133324.0, 46984.0]
[4746383236.0,9492766472.0,4746383236.0,137788.0,137788.0]
b=（-1,-1,-1,-1,-1）’;
（2）用MARLAB求低阶方程的指令A＼b得待定系数
为：
（a1,a2 ,a3,a4,a5）’
（3）a.直接法：
先将[A,b] 用高斯—约当消元法和行主元法求行最简行矩阵
R =
（a1,a2 ,a3,a4,a5）’
Matlab输入如下：
x=[5360558460628596666268894]';
y=[602611179169542349268894]';
A=[x.^22*x.*yy.^22*x2*y];
b=-ones（5,1）;
rref（[A,b]）
b.Jacobi迭代法:
取x（0）=（0,0,0,0,0）’,迭代10次，由matlab求得
（a1,a2 ,a3,a4,a5）’ =
Matlab输入如下：
x=[5360558460628596666268894]';
y=[602611179169542349268894]';
A=[x.^22*x.*yy.^22*x2*y];
b=-ones（5,1）;
X0=zeros（5,1）;
D=diag（diag（A））;
L=tril（A）-D;
U=triu（A）-D;
Bj=-D\（L+U）;
fj=D\b;
X1=A\b;
k=1;
whilek>0,
X=Bj*X0+fj;
X0=X;
if （norm（X-X1,inf）<10^-5）||（k>10）
break;
end
k=k+1;
end
c. Gauss-Seidel迭代法：
取x（0）=（0,0,0,0,0）’,迭代10次，由matlab求得
（a1,a2 ,a3,a4,a5）’
Matlab输入如下：
x=[5360558460628596666268894]';
y=[602611179169542349268894]';
A=[x.^22*x.*yy.^22*x2*y];
b=-ones（5,1）;
X0=zeros（5,1）;
D=diag（diag（A））;
L=tril（A）-D;
U=triu（A）-D;
Bg=-（D+L）\U;
fg=（D+L）\b;
X1=A\b;
k=1;
whilek>0,
X=Bg*X0+fg;
X0=X;
if （norm（X-X1,inf）<10^-5）||（k>10）
break;
end
k=k+1;
end
2.2
（1） 用Gauss列主元消去法、Gauss按比例列主元消去法、Cholesky分解求解下列线性方程组，并彼此互相验证。
（2） 判断用Jacobi迭代法、Gauss-Seidel迭代法、SOR法（分别取）解下列线性方程组的收敛性. 若收敛，再用Jacobi迭代法、Gauss-Seidel迭代法、SOR法（分别取）分别解线性方程组，并比较各种方法的收敛速度.
（3） 用Cholesky分解求解下列线性方程组
（方程组的精确解是）
解:
（1）a.Gauss列主元消去法：
步骤如下：
1 选主元
Max{a
（1）11, a
（1）21, a
（1）31, a
（1）41,}=max{1,-1,2,1}=2
即取a
（1）31作为第一次消元时的主元故需进行变换r1←→r3,原方程组变为等价方程组
执行一次消元过程，乘数（l21，l31，l41）=（-0.5，0.5，0.5）得同解方程组
2重复以上过程得
3在通过回代过程得原方程的精确解x4=2,x3=11/3,x2=1/3,x1=-8
Matlab输入如下：
A=[1-1211;-130-33;209-65;1-3-6197;];
forp=1:
3,
[Y,j]=max（abs（A（p:
4,p）））;%求主元第p列j+p-1行
C=A（p,:
）;
A（p,:
）=A（j+p-1,:
）;
A（j+p-1,:
）=C;%作行变换rp,rj+p-1;
for k=p+1:
4%消元
m=A（k,p）/A（p,p）;
A（k,p:
5）=A（k
升级会员 