数值分析上机实习报告

数值分析上机实习报告本文简介：（数值分析上机实验报告）院系：矿业学院专业：矿业工程班级：2015姓名：王学号：2015022指导教师：代第一题1.用Newton法求解方程，在（0.1，1.9）的近似根（初始近似值取为区间端点，迭代6次或误差小于0.00001）。1.1理论依据及方法应用条件Newton迭代法:由一般迭代函数，取s

数值分析上机实习报告本文内容：

（数值分析上机实验报告）

院

系：

矿业学院

专

业：

矿业工程

班

级：

2015

姓

名：

王

学

号：

2015022

指导教师：

代

第一题

1.用Newton法求解方程，在（0.1，1.9）的近似根（初始近似值取为区间端点，迭代6次或误差小于0.00001）。

1.1理论依据及方法应用条件

Newton迭代法:由一般迭代函数，取s=2时，有，可得二阶迭代序列，此种迭代法称为Newton迭代法。

定理:设函数在有限区间[a,b]上二阶导数存在，且满足条件

（Ⅰ）

（Ⅱ）在区间[a,b]上不变号;

（Ⅲ）≠0;

（Ⅳ）||/b-a≤||其中c是a,b中使min[|,]达到的一个;

则对任意时近似值x0∈[a,b],由Newton迭代过程有:

k=0,1,2…

所产生的迭代序列{x0}平方收敛于方程=0区间[a,b]上的唯一解α。

推论:设函数f(x)满足定理中条件Ⅰ,Ⅱ,Ⅲ,若选初值，使·>0,则Newton迭代过程

(k=0,1,2…)

产生的迭代序列{xk}单调收敛于=0的唯一解α。

1.2计算程序

#

include

#

include

#

include

#

include

using

namespace

std;

doublenewton

(double

a,double

b,double

eps);//牛顿迭代函数

double

newtonz

(double

x);//牛顿迭代子函数

void

main

()

{

double

a=0.1,b=1.9,eps=0.00001,*result;//初始数据

couteps)//代入a迭代计算

{

k++;

x2

=

x1;

x1

=

newtonz(x1);//调用牛顿迭代子函数

x0

=

fabs(x1-x2);

}re[0]

=

x1;

x0

=

b,k

=

0;

while

(x0>eps)//代入b迭代计算

{

k++;

x2

=

x1;

x1

=

newtonz(x1);//调用牛顿迭代子函数

x0

=

fabs(x1-x2);

}re[1]

=

x1;

return

re;

}

1.3计算结果打印

1.4

MATLAB上机程序

function

y=Newton(f,df,x0,eps,M)

d=0;

for

k=1:M

if

feval(df,x0)==0

d=2;break

else

x1=x0-feval(f,x0)/feval(df,x0);

end

e=abs(x1-x0);

x0=x1;

if

e>

x0=1.9;

>>

eps=0.00001;

>>

M=100;

>>

x=Newton(

f,df,x0,eps,M);

>>

vpa(x,7)

1.5问题讨论

1.需注意的是，要使用Newton迭代法须

满足定理中的条件Ⅰ,Ⅱ,Ⅲ，以及·>0。要用误差范围来控制循环的次数，保证循环的次数和质量，编写程序过程中要注意标点符号的使用，正确运用适当的标点符号，Newton迭代法是局部收敛的，在使用时应先确定初始值，否则所得的解可能不在所要求的范围内。

(3)因为newton法求方程是平方收敛的，所以较为精确,但是要求出函数的导数，且必须有二阶导数。

第二题

2.已知函数值如下表:

1

2

3

4

5

0

0.69314718

1.0986123

1.3862944

1.6094378

6

7

8

9

10

1.7917595

1.9459101

2.079445

2.1972246

2.3025851

=1

=0.1

试用三次样条插值求及的近似值。

2.1理论依据及方法应用条件

三次样条插值函数可定义为：对于[a,b]上的一个划分∏

a=2)

如果定义在[a,b]上的函数S(x),满足（1）.在[xi,xi+1]上为3次多项式；（2）.

S(x),S＇(x)，S＂(x)在[a,b]上连续，则称S(x)在

[a,b]上划分的3次样条函数，如果对于，还满足，，则称为的三次样条插值函数。

其基本思想是对均匀分划的插值函数的构造,三次样条函数空间中1,x,,x2,x3,(x-xj)+3为基函数,而取B样条函数Ω3(（x-xj）/h)为基函数.由于三次样条函数空间是N+3维的,故我们把分点扩大到X-1,XN+1,则任意三次样条函数可用Ω3(（x-xj）/h)线性组合来表示S(x)=

cjΩ3(（x-xj)/h)这样对不同插值问题,若能确定cj由解的唯一性就能求得S(x)。

由s(xi)=yi,I=1,2,…N

s’(x0)=y0’,s’(xN)=yN’可得

S(xi)=

cjΩ3(（xi-xj）/h)=yi

S’(x0)=1/hcjΩ3’(（x0-xj）/h)=y’0

S’(xN)=1/hcjΩ3’(（xN-xj）/h)=y’N

2.2计算程序

#include

#include

#define

N

10

/*宏定义*/

main()

{

float

s,ds,t;

float

dy0=1,dy9=0.1;

int

j;

int

x[N]={1,2,3,4,5,6,7,8,9,10};

float

y[N]={0,0.69314718,1.0986123,1.3862944,1.6094378,1.7917595,1.9459101,2.079445,2.1972246,2.3025851};

int

b[N]={2,2,2,2,2,2,2,2,2,2},h[N-1];

float

d[N],u[N-1],v[N-1],a[N-1],c[N-1],B[N],l[N-1],p[N],X[N];

for(j=1;j=0;j--)

X[j]=p[j]/B[j]-c[j]*X[j+1]/B[j];

t=4.563;

s=X[3]*pow((x[4]-t),3)/6/h[3]+X[4]*pow((t-x[3]),3)/6/h[3]+

(y[3]-X[3]*h[3]*h[3]/6)*(x[4]/h[3]-t/h[3])+

(y[4]-X[4]*h[3]*h[3]/6)*(t/h[3]-x[3]/h[3]);

/*解f(x)的值*/

ds=-X[3]*pow((x[4]-t),2)/2/h[3]+X[4]*pow((t-x[3]),2)/2/h[3]-

(y[3]-X[3]*h[3]*h[3]/6)/h[3]+(y[4]-X[4]*h[3]*h[3]/6)/h[3];

/*解f’(x)的值*/

printf(“s=%f\nds=%f\n

“,s,ds);

/*打印结果*/

}

2.3计算结果打印

2.4

MATLAB上机程序

function

Q=san(ssss,p)

Q=zeros(2,1);

x=[1;2;3;4;5;6;7;8;9;10];

y=[0;0.69314718;1.0986123;1.3862944;1.6094378;1.7917595;1.9459101;2.079445;2.1972246;2.3025851];

h=zeros(10,1);

d=zeros(10,1);

u=zeros(10,1);

v=zeros(10,1);

r=zeros(10,1);

l=zeros(10,1);

z=zeros(10,1);

m=zeros(10,1);

for

t=1:1:9;

h(t)=x(t+1)-x(t);

end

d(1)=6/h(1)*((y(2)-y(1))/h(1)-1);

d(10)=6/h(9)*(0.1-(y(10)-y(9))/h(9));

for

t=1:1:8

u(t+1)=h(t)/(h(t)+h(t+1));

v(t+1)=1-u(t+1);

d(t+1)=6/(h(t)+h(t+1))*((y(t+2)-y(t+1))/(x(t+2)-x(t+1))-(y(t+1)-y(t))/(x(t+1)-x(t)));

end

u(10)=1;v(1)=1;r(1)=d(1);

for

t=2:1:10

l(t)=u(t)/r(t-1);

r(t)=d(t)-l(t)*v(t-1);

end

z(1)=d(1);

for

t=2:1:10

z(t)=d(t)-l(t)*z(t-1);

end

m(10)=z(10)/r(10);

for

t=9:-1:1

m(t)=(z(t)-v(t)*m(t+1))/r(t);

end

for

t=1:1:10

if

p>=t

y=pow(3,x)*pow(x,1.4)*(5*x+7)*sin(x*x);

return(y);

}

main()

{

float

T[N+1][N+1],h[N+1],a=1,b=3,m[N+1];

int

i,l;

T[1][0]=(b-a)*(f(a)+f(b))/2;

l=1;

while(leps)

J=J+1;h=h/2;S=0;

for

i=1:n

x=a+h*(2*i-1);

S=S+feval(f,x);

end

R(J+1,1)=R(J,1)/2+h*S;

for

k=1:J

R(J+1,k+1)=(4^k*R(J+1,k)-R(J,k))/(4^k-1);

end

err=abs(R(J+1,J+1)-R(J+1,J));

n=2*n;

end

R;

T=R(J+1,J+1);

format

long

[emailprotected](x)(3.^x)*(x.^1.4)*(5*x+7)*sin(x*x);

[T,n]=mromb(f,1,3,1.e-5)

3.5问题讨论

1、Romberge算法的优点是:

a、把积分化为代数运算,而实际上只需求T1(i),以后用递推可得。

b、算法简单且收敛速度快,一般4或5次即能达到要求。

c、节省存储量，算出的可存入。

2、Romberge算法的缺点是:

a、对函数的光滑性要求较高。

b、计算新分点的值时,这些数值的个数成倍增加。

第四题

4.用定步长四阶Runge-Kutta法求解,打印,,,,4.1理论依据及方法应用条件

Runge-Kutta法的基本思想:不是按Taylor公式展开，而是先写成处附近的值的线性组合（有待定系数）,再按Taylor公式展开，然后确定待定常数。

四阶古典Runge-Kutta公式:

4.2计算程序

#include

int

main()

{

int

i;

double

h=0.0005;

double

k1,k2,k3,k4;

double

y1=0.0,y2=0.0,y3=0.0;

for(i=1;i0,这样只需要把上式右端第二项取为2sign(a21)a21

s1即可，从而

归纳起来算法步骤为：

（2）程序主体

#include

#include

void

main(

)

{

int

i,j,m,r,sign;

double

A0[9][9],s,z,p,n,v,h,l,y[9],u[9],k,q[9],X[9],x[9]={0,0,0,0,0,0,0,0,0},B[9][9],g[9];

double

A[9][9]=

{

{12.38412,2.115237,-1.061074,1.112336,-0.113584,0.718719,1.742382,3.067813,-2.031743},{2.115237,19.141823,-3.125432,-1.012345,2.189736,1.563849,-0.784165,1.112348,3.123124},{-1.061074,-3.125432,15.567914,3.123848,2.031454,1.836742,-1.056781,0.336993,-1.010103},{1.112336,-1.012345,3.123848,27.108437,4.101011,-3.741856,2.101023,-0.71828,-0.037585},{-0.113584,2.189736,2.031454,4.101011,19.897918,0.431637,-3.111223,2.121314,1.784317},{0.718719,1.563849,1.836742,-3.741865,0.431637,9.789365,-0.103458,-1.103456,0.238417},{1.742382,-0.784165,-1.056781,2.101023,-3.111223,-0.103458,14.713846,3.123789,-2.213474},{3.067813,1.112348,0.336993,-0.71828,2.121314,-1.103456,3.123789,30.719334,4.446782},{-2.031743,3.123124,-1.010103,-0.037585,1.784317,0.238417,-2.213474,4.446782,40.00001}

};

double

b[9]={2.1874369,33.992318,-25.173417,0.84671695,1.784317,-86.612343,1.1101230,4.719345,-5.6784392};

for(i=0;i0)

sign=1;

else

if(A[r+1][r]k

A([k

p],:)=A([p

k],:);

b([k

p],:)=b([p

k],:);

end

m=A(k+1:n,k)/A(k,k);

A(k+1:n,k+1:n)=A(k+1:n,k+1:n)-m*A(k,k+1:n);

b(k+1:n)=b(k+1:n)-m*b(k);

A(k+1:n,k)=zeros(n-k,1);

if

flag~=0,Ab=[A,b],end

end

x=zeros(n,1);

x(n)=b(n)/A(n,n);

for

k=n-1:-1:1

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

end

format

long

A=[12.38412,2.115237,-1.061074,1.112336,-0.113584,0.718719,1.742382,3.067813,-2.031743;

2.115237,19.141823,-3.125432,-1.012345,2.189736,1.563849,-0.784165,1.112348,3.123124;

-1.061074,-3.125432,15.567914,3.123848,2.031454,1.836742,-1.056781,0.336993,-1.010103;

1.112336,-1.012345,3.123848,27.108437,4.101011,-3.741856,2.101023,-0.71828,-0.037585;

-0.113584,2.189736,2.031454,4.101011,19.897918,0.431637,-3.111223,2.121314,1.784317;

0.718719,1.563849,1.836742,-3.741856,0.431637,9.789365,-0.103458,-1.103456,0.238417;

1.742382,-0.784165,-1.056781,2.101023,-3.111223,-0.103458,14.713846,3.123789,-2.213474;

3.067813,1.112348,0.336993,-0.71828,2.121314,-1.103456,3.123789,30.719334,4.446782;

-2.031743,3.123124,-1.010103,-0.037585,1.784317,0.238417,-2.213474,4.446782,40.00001];

b=[2.1874369,33.992318,-25.173417,0.84671695,1.784317,-86.612343,1.1101230,4.719345,-5.6784392]

;

x=mgauss2(A,b);x=x

（5）问题讨论

1、矩阵A是严格对角占优阵，而且还是实数对称阵，满足Householder变换的条件，并且得到的B就是与A相似的三对角阵。

2、经过Householder变换矩阵有良好的性质，便于计算分析。