贵州大学2014级数值分析上机实验

  1. 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
  2. 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
  3. 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。

一、用Newton法求方程:x7-28x4+14=0 在(0.1,1.9)中的近似根(初始近似值取为区间

端点,迭代6次或误差小于0.00001.)

程序编写:

#include

#include

main()

{

double x0,x1,y,y0,y1;

printf("x0=");

scanf("%lf",&x0);

do

{

y0=pow(x0,7)-28*pow(x0,4)+14;

y1=7*pow(x0,6)-112*pow(x0,3);

x1=x0-y0/y1;

y=fabs(x1-x0);

x0=x1;

}

while(y>=0.00001);

printf("x=%lf\n",x1);

}

二、已知函数值如下表:

#include

#include

#define N 9

void main()

{

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

y[N+1]={0,0.69314718,1.0986123,1.3862944,1.6094378,1.7917595,1.9459101,2.079445,2.1 972246,2.3025851},

h[N+1],d[N+1],a[N+1],c[N+1],b[N+1]={2,2,2,2,2,2,2,2,2,2},s[N+1],t[N+1],l[N+1],M[N+1],f,f1;

int i;

for(i=1;i<=N;i++)

h[i-1]=x[i]-x[i-1];

d[0]=6/h[0]*((y[1]-y[0])/h[0]-1);

d[N]=6/h[N-1]*(0.1-(y[N]-y[N-1])/h[N-1]);

for(i=1;i<=N-1;i++)

{

d[i]=6/(h[i-1]+h[i])*((y[i+1]-y[i])/h[i]-(y[i]-y[i-1])/h[i-1]);

a[i]=h[i-1]/(h[i-1]+h[i]);

c[i]=1-a[i];

}

c[0]=1;a[N]=1;s[0]=b[0];t[0]=d[0];

for(i=0;i<=N-1;i++)

{

l[i+1]=a[i+1]/s[i];

s[i+1]=b[i+1]-l[i+1]*c[i];

t[i+1]=d[i+1]-l[i+1]*t[i];

}

M[N]=t[N]/s[N];

for(i=N-1;i>=0;i--)

M[i]=(t[i]-c[i]*M[i+1])/s[i];

f=M[3]*pow((x[4]-4.563),3)/6/h[3]

+M[4]*pow((4.563-x[3]),3)/6/h[3]

+(y[3]-M[3]*h[3]*h[3]/6)*(x[4]-4.563)/h[3]

+(y[4]-M[4]*h[3]*h[3]/6)*(4.563-x[3])/h[3];

f1=-3*M[3]*pow((x[4]-4.563),2)/6/h[3]

+3*M[4]*pow((4.563-x[3]),2)/6/h[3]

-(y[3]-M[3]*h[3]*h[3]/6)/h[3]

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

printf("f(4.563)=%lf f'(4.563)=%lf\n",f,f1); }

三,用Romberg算法求

3 1.42

1

3(57)sin

x x x x dx

+

⎰(允许误差ε=0.000 01)。

程序:

#include

#include

float f(float x)

{

float f=0.0;

f=pow(3.0,x)*pow(x,1.4)*(5*x+7)*sin(x*x);

return (f);

}

main()

{

int i=1,j,k,n=12;

float T[12],a=1.0,b=3.0,s=0.0; T[0]=0.5*(b-a)*(f(a)+f(b)); for(j=1;j

s+=f(a+(2*k-1)*(b-a)/pow(2,j));

T[j]=0.5*(T[j-1]+(b-a)*s/pow(2,j-1)); s=0.0;

}

T[11]=(4*T[1]-T[0])/(float)3;

for(;fabs(T[11]-T[0])>0.00001;i++) { T[0]=T[11];

for(j=1;j

printf("%f\n",T[11]);

}

4. 用定步长四阶Runge-Kutta 求解

⎪⎪⎪⎪⎩⎪⎪⎪⎪⎨

⎧===--===0

)0(0)0(0)0(10010001000//1/321

3

2332

1y y y y y dt dy y

dt dy dt dy h =0.0005,打印y i (0.025) , y i (0.045) , y i (0.085) , y i (0.1) ,(i =1,2,3)

相关文档
最新文档