数值分析课程设计

bd336x280();=22-jpg_6_0_______-50-0-284-50.jpg" alt="数值分析课程设计" />

学院：数学学院

姓名：王敏

班级：数学与应用数学一班 学号：20096269 时间：2011/6/18

实验一（2题)

1.1 水手、猴子和椰子问题：五个水手带了一只猴子来到南太平洋的一个荒岛

上，发现那里有一大堆椰子。由于旅途的颠簸，大家都很疲惫，很快就入睡了。第一个水手醒来后，把椰子平分成五堆，将多余的一只给了猴子，他私藏了一堆后便又去睡了。第二、第三、第四、第五个水手也陆续起来，和第一个水手一样，把椰子分成五堆，恰多一只猴子，私藏一堆，再去入睡，天亮以后，大家把余下的椰子重新等分成五堆，每人分一堆，正好余一只再给猴子，试问原先共有几只椰子？

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

解：首先分析椰子数目的变化规律，设最初的椰子数为p

后剩余的椰子数目，则根据问题有

，即第一个水手所处

理之前的椰子数，用p 1、p 2、p 3、p4、p 5 分别表示五个水手对椰子动了手脚以

pk 1

45

(pk 1)(k 0,1,2,3,4)

再用x表示最后每个水手平分得到的椰子数，于是有

x

15

(p5 1)

所以

p5 5x 1

利用逆向递推的方法，有

pk

54

pk 1 1,(k 4,3,2,1,0)

有了逆向递推关系式，求解这一问题似乎很简单，但由于椰子数为一正整数，用任意的x作为初值递推出的p0数据不一定是合适的。 这里用 for 循环语句结合 break 语句来寻找合适的 x 和 p0 ，对任意的 x 递推计算出 p0 ，当计算结果为正整数时，结果正确，否则选取另外的 x 再次重新递推计算，直到计算出的结果 p0 为正整数为止。

程序如下 ：

n=input('input ：n') for x=1:n p=5*x+1; for k=1:5 p=5*p/4+1; end

if p==fix(p), break, end end

disp([x,p])

运行这段程序后，屏幕出现要求从键盘输入 x 数据的信息input1023后，MATLAB计算出合适的 x 和 p0 的值。 运行结果如下： input ：n 1023 n =

1023

1023 15621

1.2n

设，I

1x

n

5 x

dx

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

I1n 5In 1

n

(n 1,2, 20)

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

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

I1n 1

5In

15n

(n 30,29, ,3,2)

计算从I1到I20的近似值；

（3）分析所得结果的可靠性以及出现这种现象的原因。

解：(1)递推公式程序如下：

function dtgs(n) syms x i I; f=1/(5+x); I=int(f,x,0,1);

n，输入

for i=1:n I=-5*I+1/i; vpa(I) end end

在MATLAB命令窗口输入：dtgs(20)

运行结果为： ans =

0.088392216030226868941409874227427 ans =

0.058038919848865655292950628862866 ans =

0.043138734089005056868580189019004 ans =

0.034306329554974715657099054904979 ans =

0.028468352225126421714504725475105 ans =

0.024324905541034558094143039291144 ans =

0.021232615151970066672141946401424 ans =

0.018836924240149666639290267992878 ans =

0.016926489910362777914659771146719 ans =

0.015367550448186110426701144266407 ans =

0.014071338668160356957403369577056 ans =

0.012976639992531548546316485448045 ans =

0.012039876960419180345340649682855 ans =

0.011229186626475526844725323014428 ans =

0.010520733534289032443040051594081 ans =

0.0098963323285548377847997420264414 ans =

0.0093418677689905169583542310372102 ans =

0.00884621671060297076378440027414 ans =

0.0084004953943535672337095771302562 ans =

0.0079975230282321638314521117637798 （2）程序为：

function dts2(n) syms x i I; f=x^30/(5+x); I=int(f,x,0,1); for i=n:-1:1 I=-1/5*I+1/(5*i); vpa(I) end end

在MATLAB命令窗口输入：dtgs2(30) 运行的结果为

ans =

0.0055857400736385308149005088362848 ans =

0.0057794037094102248716779050494163 ans =

0.0059869764009750978827951756622988 ans =

0.0062100121272123878308468036954157 ans =

0.0064503052668652147415239650216708 ans =

0.0067099389466269570516947040698534 ans =

0.0069913455440079419229941995105531 ans =

0.0072973830651114550936620341586536 ans =

0.0076314324778867998903584895694758 ans =

0.0079975230282321638314521085326055 ans =

0.0084004953943535672337095779380497 ans =

0.0088462167106029707637844000721916 ans =

0.0093418677689905169583542310876973 ans =

0.0098963323285548377847997420201305

ans =

0.010520733534289032443040051595658 ans =

0.011229186626475526844725323014033 ans =

0.012039876960419180345340649682905 ans =

0.012976639992531548546316485448033 ans =

0.014071338668160356957403369577059 ans =

0.015367550448186110426701144266406 ans =

0.016926489910362777914659771146719 ans =

0.018836924240149666639290267992878 ans =

0.021232615151970066672141946401424 ans =

0.024324905541034558094143039291144 ans =

0.028468352225126421714504725475105 ans =

0.034306329554974715657099054904979 ans =

0.043138734089005056868580189019004 ans =

0.058038919848865655292950628862866 ans =

0.088392216030226868941409874227427 ans =

0.18232155679395462621171802515451

（3）第一小题：可以看到虽然初始误差很小，但是上述算法的误差是逐步扩大的，算法是不稳定的，因此结果不可靠。产生错误的原因：递推公式每次将前一步的误差扩大5倍，这会使误差迅速扩大，而且舍入误差还会影响结果，所以不稳定。第二小题：有数据可以看出，该算法得出的的精度很高，因为误差的每一步都是缩小的，所以算法稳定。(该例子告诉我们，用数值方法解决问题时一定要选择数值稳定的算法)

实验二（3题）

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

20 A 1

2

28 3

3

1 15

消元过程可将矩阵A化为上三角矩阵U，试求出消元过程所用的乘数l21、l31、l32并以如下格式构造下三角矩阵L和上三角矩阵U

1

L l21

l31

1l32

20

,U

1

2u22

(1)

(1) u23 (2)u33

3

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

解：将矩阵A进行直接LU分解的程序

function [L,U]=myLU(A) [n n]=size(A); L=zeros(n,n); U=zeros(n,n); for i=1:n L(i,i)=1; end for k=1:n for j=k:n

U(k,j)=A(k,j)-sum(L(k,1:k-1).*U(1:k-1,j)'); end for i=k+1:n

L(i,k)=(A(i,k)-sum(L(i,1:k-1).*U(1:k-1,k)'))/U(k,k); end end L(2,1) L(3,1) L(3,2) end

在MATLAB命令窗口输入：

A=[20 2 3;1 8 1;2 -13 15];[L,U]=myLU(A)

运行的结果为

ans = 0.0500 ans = 0.1000

ans = -1.6709

L =

1.0000 0 0 0.0500 1.0000 0 0.1000 -1.6709 1.0000 U =

20.0000 2.0000 3.0000 0 7.9000 0.8500 0 0 16.1203

经验证可得A=LU,且得：

l32 1.6709

l21 0.0500 l31 0.1000

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

1 0 0 b1 0,b2 1,b3 0

0 0 1

分别求解方程组

AX b1,AX b2,AX b3

由于三个方程组系数矩阵相同，可以将分解后的矩阵重复使用。对第一个方程组，由于A=LU，所以先求解下三角方程组LY b1,再求解上三角方程组

UX Y

，则可得逆矩阵的第一列列向量；类似可解第二、第三方程组，得逆矩

阵的第二列列向量的第三列列向量。由三个列向量拼装可得逆矩阵A 1。

解：矩阵分解程序同上一题，下面解方程AX=b1：

在MATLAB命令窗口输入：

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

b=[1 0 0]'; [L,U]=lu(A) X=U\(L\b) 运行的结果为

X= 0.0517

-0.0055 -0.0080

同理可以求得AX=b2，AX=b3的解：

X= X=

-0.0164 -0.0093

0.1237 -0.0072

0.0269 0.0665 所以A的逆矩阵为

-0.0164 0.0517

-0.0055 0.1237 0.0269 -0.0080

-0.00

-0.0072 0 .06 65

2.3验证希尔伯特矩阵的病态性：对于三阶矩阵

1

H 1/2

1/3

1/21/31/4

1/3

1/4 1/5

取右端向量b [11/613/12（1）向量X [x1

x2

T

47/60]1

1]

T

，验证：

是方程组HX b的准确解；

x3] [1

T

（2）取右端向量b的三位有效数字得b [1.831.080.783]T，求方程组的准确解X ，并与X的数据[111]T作比较 。说明矩阵的病态性。

解：希尔伯特矩阵的病态性程序：

(1)function [RA,RB,n,x]=btjz(A,b)

B=[A b]; n=length(b); RA=rank(A); RB=rank(B);

zhica=RB-RA; if zhica>0 return end

if RA==RB if RA==n

x=zeros(n,1); C=zeros(1,n+1); for p=1:n-1

[Y,j]=max(abs(B(p:n,p))); C=B(p,:); B(p,:)=B(j+p-1,:);B(j+p-1,:)=C; for k=p+1:n

m=B(k,p)/B(p,p);

B(k,p:n+1)=B(k,p:n+1)-m*B(p,p:n+1); end end

b=B(1:n,n+1);A=B(1:n,1:n);x(n)=b(n)/A(n,n); for q=n-1:-1:1

x(q)=(b(q)-sum(A(q,q+1:n)*x(q+1:n)))/A(q,q);

end else end end

在MATLAB命令窗口输入：

A=[1 1/2 1/3;1/2 1/3 1/4;1/3 1/4 1/5];

b=[11/6;13/12;47/60];[RA,RB,n,x]=btjz(A,b) 运行结果： RA = 3 RB = 3 n = 3 x =

1.0000 1.0000 1.0000

（2）在MATLAB命令窗口输入：

A=[1 1/2 1/3;1/2 1/3 1/4;1/3 1/4 1/5]; b=[1.83;1.08;0.783];[RA,RB,n,x]=btjz(A,b) 运行结果：

RA = 3 RB = 3 n = 3 x =

1.0800

0.5400 1.4400

由于detA的值相对较小，所以是病态的。

实验三 （2题）

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

y0(x) sinx,x [0, ]y1(x) x,x [0, /2]y2(x) x x/6,x [0, /2]

y3(x) x x/6 x/120,x

[0, /2]

3

5

3

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

解：画图程序：（在MATLAB窗口输入：）

1、 x=[0:0.5:360]*pi/180; plot(x,sin(x)) ;

2、 x=[0:0.5:180]*pi/180; plot(x,x) ;

3、 syms x ;

ezplot('x-x^2/6');

4、 syms x ;

ezplot('x-x^2/6+x^5/120');

运行结果：

3.2绘制飞机的降落曲线

一架飞机飞临北京国际机场上空时，其水平速度为540km/h，飞行高度为1 000m。飞机从距机场指挥塔的横向距离12 000m处开始降落。根据经验，一架水平飞行的飞机其降落曲线是一条三次曲线。建立直角坐标系，设飞机着陆点为原点O，降落的飞机为动点P(x,y)，则x表示飞机距指挥塔的距离，y表示飞机的飞行高度，降落曲线为

y(x) a0 a1x a2x a3x

2

3

该函数满足条件：

y(0) 0,y(12000) 1000y'(0) 0,y'(12000) 0

（1）试利用y(x)满足的条件确定三次多项式中的四个系数； （2）用所求出的三次多项式函数绘制出飞机降落曲线。

解:（1）.程序如下：

function f=Hermite(x,y,y_1,x0) syms t; f=0.0; n=length(x) for i=1:n; h=1.0; a=0.0; for j=1:n if j~=i

h=h*(t-x(j))^2/((x(i)-x(j))^2); a=a+1/(x(i)-x(j)); end end

f=f+h*((x(i)-t)*(2*a*y(i)-y_1(i))+y(i)); if i==n

if nargin==4 f=sub(f,'t',x0); else

f=vpa(f,6); end end end

在MATLAB命令窗口输入：

x=[0 1200];y=[0 1000];y_1=[0 0];f=Hermite(x,y,y_1) 运行结果为： n = 2

f =-6.94444*10^(-7)*t^2*(1.66667*t - 3000.0) 由函数f即可得系数.

（2）画图程序为（在MATLAB窗口输入：） syms x ;

ezplot('-6.94444*10^(-7)*t^2*(1.66667*t - 3000.0)'); 图像：

实验四（1题）

4.1

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

现象：集成电路上可容纳的零件数量，每隔一年半左右就会增长一倍，性能也提升一倍。因而发表论文，提出了大大有名的摩尔定律（Moore’s Law），并预测未来这种增长仍会延续下去。下面数据中，第二行数据为晶片上晶体数目在不同年代与1959年时数目比较的倍数。这些数据是推出摩尔定律的依据：

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

解：利用最小二乘法，其程序为：

function zxec(x,y,n) y=y';

for i=1:length(x) for j=1:n+1

A(i,j)=x(i)^(j-1); end end

[L,U]=lu(A'*A); alpha=U\(L\(A'*y))

在MATLAB命令窗口输入：

x=[1,4,5,6,7];y=[1,3,4,5,6];zxec(x,y,1)

运行结果:

alpha = -0.0189

0.8302

在MATLAB命令窗口输入：

x=[1,4,5,6,7];y=[1,3,4,5,6];zxec(x,y,2) alpha = 0.4961 0.4437

0.0497 函数表达式：

f=0.496075353218219x^2+0.443746729461008x+0.049712192569336 画图的程序： syms x ;

ezplot('0.496075353218219*x^2+0.443746729461008*x+0.049712192569336')

图像:

实验五（3题）

10

5.1用几种不同的方法求积分

41 x

2

dx

的值。

（1）牛顿-莱布尼茨公式；（2）梯形公式；（3）辛卜生公式；（4）复合梯形公式。 解：

（1）由牛顿-莱布尼茨公式公式解得

10

41+x

0= -0=3.1415926 2

1

（2）梯形公式的程序如下(在MATLAB命令窗口输入)：

y=0;

for x=0:1

y=y+4/(1+x^2);

end (1/2)*y 运行结果如下:

ans =

3

(3) 辛卜生公式的程序如下

function fhSimpson(a,b,n)

h=(b-a)/n;

x=a:h:b;

S1=0; S2=0; for i=1:n

S1=S1+f(a+(i-1/2)*h); end

for i=1:n-1 S2=S2+f(a+i*h); end

format long

S=(b-a)/(6*n)*(f(a)+4*S1+2*S2+f(b))

现新建f.m文件： function r=f(x) r=4/(1+x^2);

在MATLAB命令窗口输入：

fhSimpson(0,1,3)

运行结果如下: S =

3.141591780936043

（4）复合梯形公式的程序如下 function fhtx(a,b,n) h=(b-a)/n;

D=f(a)/2+f(b)/2; for i=1:n-1

D=D+f(a+i*h); end

Tn=h*D

现新建f.m文件： function r=f(x) r=4/(1+x^2);

在MATLAB命令窗口输入：

fhtx(0,1,8) 运行结果如下: Tn =

3.138988494491089

5.2

设X为标准正态随机变量，即X~N（0，1）。现分别取

30个不同的概率值；P X ua ，并

u 0.1,0.2,0.3, ,2.9,3.0，试设计算法计算

将计算结果与概率论教科书中的标准正态分布函数表作比较。

（提示： P

X ua

ua

e

x

2

2

dx

4

ua

e

x

2

2

dx)

解:程序如下(在MATLAB命令窗口中输入：)

for u=0.1:0.1:3.0

p=int(1/sqrt(2*pi)*exp(-x^2/2),u,4); vpa(p) end

运行结果如下： ans =

0.46014049148113792735987262882212 ans =

0.4207086193190638833193970255296 ans =

0.3820569065692142666406029409212 ans =

0.34454658714784273486685485083196 ans =

0.30850586748415379571444540073965 ans =

0.27422144650824047750464919332639 ans =

0.24193198098123990994240483988631 ans =

0.21182372734156357888762186243779 ans =

0.18402845410492638012982327033221 ans =

0.15862358268962394140326477999267 ans =

0.13563438970454956372534181236522 ans =

0.11503799897987515528777887093705 ans =

0.09676881334377721927623074191764 ans =

0.080724987991937931618097576885656 ans =

0.066775530027024950254933222030685 ans =

0.054767620457724877460738469926159 ans =

0.044533791516709922348355640672936 ans =

0.035898647871092686281781901339917 ans =

0.028684888574168681272125090189663 ans =

0.022718460706346088698327857615476 ans =

0.017832749320983437976742564767419 ans =

0.013871776271665491560116833525591 ans =

0.01069243877984268613909448646859 ans =

0.0081658646827630100332823234341801 ans =

0.006177994083943015631684495659137 ans =

0.0046295167818856306189611919583325 ans =

0.0034353025612075487893030127714414 ans =

0.0025234590885948130379266044246572 ans =

0.0018341420585509181436415648772698 ans =

0.0013182267897969746877521249991372

经比较知计算结果与概率论教科书中的标准正态分布函数表的知近似。

x cost

5.5 求空间曲线L： y sint

z 2 cost sint

弧长公式为

L

2

解：先计算可知积分函数为(2-sin(2*t))^(1/2)，

程序如下(在MATLAB命令窗口中输入：)

syms t s;

s=int((2-sin(2*t))^(1/2),t,0,2*pi); vpa(s)

运行结果如下： ans =

8.737752570984804741619351957708

实验六（1题）

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

(1)y'

0.9y1 2x

xy1 x

2

,y(0) 1;x [0,1]

,y(0) 2;x [0,1]

(2)y'

.

解: 欧拉公式程序如下:

function euler(a,b,h,alpha) x=a:h:b; n=(b-a)/h; y(1)=alpha; for i=2:n+1

y(i)=y(i-1)+h*f(x(i-1),y(i-1));

sprintf('x=%3.1f,y=%9.7f',x(i),y(i)) end

r =

ans =

x=0.4,y=0.0000000 r =

ans =

x=0.5,y=0.0000000 r=

ans =

x=0.6,y=0.0000000 r =

ans =

x=0.7,y=0.0000000 r =

ans =

x=0.8,y=0.0000000 r =

ans =

x=0.9,y=0.0000000 r =

ans =

x=1.0,y=0.0000000

先建立名为f.m的文件 function r=f(x,y) r=(0.9*y)/(1+2*x)

在MATLAB命令窗口输入：

euler(0,1,0.1,0)

运行结果 r =

ans =

x=0.1,y=0.0000000 r =

ans =

x=0.2,y=0.0000000 r =

ans =

x=0.3,y=0.0000000

四阶龙格-库塔法程序如下： function rk4(a,b,h,alpha) t=a:h:b; n=(b-a)/h; y(1)=alpha; for i=2:n+1

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);

sprintf('t=%3.1f,y=%9.7f',t(i),y(i)) end

先建立名为f.m的文件 function r=f(x,y) r=-(t*y)/(1+t^2)

在MATLAB命令窗口输入： rk4(0,1,0.1,0) 运行结果 r =

0 r = 0 r =

0 r = 0 ans =

t=0.1,y=0.0000000 r =

0 r =

0 r =

r = 0 ans =

t=0.2,y=0.0000000 r = 0 r =

0 r =

0 r = 0 ans =

t=0.3,y=0.0000000 r =

0 r =

0 r =

0 r =

0 ans =

t=0.4,y=0.0000000 r = 0 r = 0 r = 0 r = 0 ans =

t=0.5,y=0.0000000 r = 0 r = 0 r = 0 r = 0 ans =

t=0.6,y=0.0000000

r = 0 r = 0 r = 0 r = 0 ans =

t=0.7,y=0.0000000 r = 0 r = 0 r = 0 r = 0 ans =

t=0.8,y=0.0000000

实验七(2题)

r =

0 r =

0 r = 0 r =

0 ans =

t=0.9,y=0.0000000 r =

0 r =

0 r =

0 r = 0 ans =

t=1.0,y=0.0000000

7.1

用二分法求下列方程在指定区间[a，b]上的实根近似值：（要求误差不超过 (2)xsin x=1，[a,b]=[0,1.5].

0.01） (1)x-sinx-1=0，[a,b]=[1,3]；

解：二分法程序如下：

function [x,k]=eff(fun,a,b,eps) if nargin<4 eps=le-5; end

fa=feval(fun,a); fb=feval(fun,b); if fa*fb>0 return; end k=0;

while abs(b-a)/2>eps

x=(a+b)/2; fx=feval(fun,x); if fx*fa<0 b=x; fb=fx; else

a=x;fa=fx;

本文来源：https://www.bwwdw.com/article/qt01.html