最新版数值计算实验报告Word文档格式.docx

最新版数值计算实验报告Word文档格式.docx

《最新版数值计算实验报告Word文档格式.docx》由会员分享，可在线阅读，更多相关《最新版数值计算实验报告Word文档格式.docx（27页珍藏版）》请在冰豆网上搜索。

姓名：

宋元台学号：

实验成绩：

1、实验目的及要求

实验目的：

掌握Newton插值多项式的算法，理解Newton插值多项式构造过程中基函数的继承特点，掌握差商表的计算特点。

实验要求：

1.给出Newton插值算法

2.用C语言实现算法

二、实验内容

三、实验步骤（该部分不够填写.请填写附页）

1.算法分析：

下面用伪码描述Newton插值多项式的算法：

Step1输入插值节点数n，插值点序列{x（i）,f（i）},i=1,2,……,n,要计算的插值点x.

Step2形成差商表

fori=0ton

forj=ntoif（j）=（（f（j）-f（j-1）/（x（j）-x（j-1-i））;

Step3置初始值temp=1,newton=f（0）

Step4fori=1ton

temp=（x-x（i-1））*temp/*由temp（k）=（x-x（k-1））*temp（k-1）形成（x-x（0）.....（x-x（i-1）*/

Newton=newton+temp*f（i）;

Step5输出f（x）的近似数值newton（x）=newton.

2.用C语言实现算法的程序代码

#include<

stdio.h>

#defineMAX_N20

typedefstructtagPOINT

{

doublex;

doubley;

}POINT;

intmain（）

intn;

inti,j;

POINTpoints[MAX_N+1];

doublediff[MAX_N+1];

doublex,tmp,newton=0;

printf（"

\nInputnvalue:

"

）;

scanf（"

%d"

&

n）;

if（n>

MAX_N）

{

theinputnislargerthanMAX_N,pleaseredefinetheMAX_N.\n"

return1;

}

if（n<

=0）

pleaseinputanumberbetween1and%d.\n"

MAX_N）;

nowinputthe（x_i,y_i）i=0,...%d\n"

n）;

for（i=0;

i<

=n;

i++）

pleaseinputx（%d）y（%d）\n"

i,i）;

%lf%lf"

points[i].x,&

points[i].y）;

nowinputthevalueofx:

\n"

%lf"

x）;

i++）diff[i]=points[i].y;

n;

for（j=n;

j>

i;

j--）

{

diff[j]=（diff[j]-diff[j-1]）/（points[j].x-points[j-1-i].x）;

}

tmp=1;

newton=diff[0];

tmp=tmp*（x-points[i].x）;

newton=newton+tmp*diff[i+1];

newton（%f）=%f\n"

x,newton）;

return0;

}

程序运行结果过程及截图为：

2014至2015第1学期2014年12月1日

2实验项目Romberg算法指导教师：

一、实验目的及要求

熟练掌握数值积分的基本思想和动手推导积分公式的能力，并写出自己设计的变步长算法，用此算法解决实际问题。

主要锻炼分析问题及对数值方法产生的误差在实际执行算法中的变化，体会事前误差和事后误差的不同之处。

写出Romberg算法的描述，并用此算法计算下列定积分的值，要求满足给定误差。

二、实验内容。

（1）用变步长的Romberg方法计算定积分：

给定误差限为

下面用伪码描述Romberg的算法描述：

Sstep1输入区间端点a,b,精度控制值e，循环次数M，定义函数f（x），取n=1，

h=b-a；

Step2R1,1=（f（a）+f（b））*h/2;

Step3fork=2toM

{Rk,1=（R（k-1）,1+h（k-1）*/**/

forj=2tok

if退出循环

Step4输出.

math.h>

#definef（x）（sin（x））

#defineN_H20

#defineMAXREPT10

#definea1.0

#defineb2.0

#defineepsilon0.00001

doublecomputeT（doubleaa,doublebb,longintn）

inti;

doublesum=0.0;

doubleh=（bb-aa）/n;

for（i=1;

sum+=f（aa+i*h）;

sum+=（f（aa）+f（bb））/2;

return（h*sum）;

voidmain（）

longintn=N_H,m=0;

doubleT[MAXREPT+1][2];

T[0][1]=computeT（a,b,n）;

n*=2;

for（m=1;

m<

MAXREPT;

m++）

for（i=0;

m;

T[i][0]=T[i][1];

for（i=1;

=m;

T[i][1]=T[i-1][1]+（T[i-1][1]-T[i-1][0]）/（pow（2,2*m）-1）;

if（（T[m-1][1]<

T[m][1]+epsilon）&

&

（T[m-1][1]>

T[m][1]-epsilon））

theintegrateis%lf\n"

T[m][1]）;

return;

returnnosolved...\n"

}_

3.运行结果

数值计算方法与算法专业:

3实验项目：

求拟合算法指导教师：

实验成绩：

学习了曲线拟合的最小二乘法，会求拟合曲线，线性拟合和二次拟合函数，还有形如aexp（bx）的曲线拟合，实验主要是求形如aexp（bx）的曲线拟合，并理解体会其意义。

1.形如aexp（bx）曲线拟合的算法；

用形如p（x）=aexp（bx）的函数拟合给定（）

拟合函数aexp（bx）的算法描述：

Step1输入m值，及（xi,yi）,i=1,2,...,m;

Step2解方程组

m

a=

b

Step3输出p（x）=aexp（bx）即所求拟合函数

#defineMAX_N25

intm;

POINTpoints[MAX_N];

staticdoubleu11,u12,u21,u22,c1,c2;

doubleA,B,tmp;

\nInputmvalue"

m）;

if（m>

=MAX_N）

TheinputmislargerthanMAX_N,pleaseredefinetheMAX_N.\n"

if（m<

Pleaseinputanumberbetween1and%d"

Nowinputthe（x_i,y_i）,i=1,...,%d:

m）;

tmp）;

points[i].x=tmp;

points[i].y=tmp;

/*inputx_i,y_i*/

u21+=points[i].x;

u22+=points[i].x*points[i].x;

c1+=log（points[i].y）;

c2+=points[i].x*log（points[i].y）;

u12=u21;

u11=m;

A=（c1*u22-c2*u12）/（u11*u22-u12*u21）;

B=（c1*u21-c2*u11）/（u21*u12-u22*u11）;

Solve:

p（x）=%lfexp（%lfx）\n"

exp（A）,B）;

3.实验输出结果：

2014至2015第1学期2014年12月1日

数值计算方法与算法专业:

4实验项目弦截断法指导教师：

学习什么是弦截法，了解弦截法的用法，并学会用弦截法求方程的根。

用弦截法求f（x）在x0,x1附近的根

弦截法的算法描述：

Sstep1定义函数f（x0,输入控制精度epsilon，迭代初始值x_k1,x_k

计算f1:

=f（x_k1）!

x_k1,x_k2表示x（k-1）,xk

Step2fork=2toMAXREPT

2.1f2:

=f（x_k2）

2.1x_k:

=x_k2-f2（x_k2-x_k1）/（f2-f1）!

x_k表示x（k+1）

2.3if（|x_k-x_k2|<

epsilon）OR（|f（x_k）|<

epsilon）

Then{输出满足给定精度的近似解x_k,结束}

2.4f1：

=f2!

为下一次迭代准备数值

x_k1:

=x_k2

X_k2:

=x_k

Step3输出：

在初始值x_k1,x_k2附近f（x）无根

4.用C语言实现算法的程序代码

#definef（x）（x*x*x-7.7*x*x+19.2*x-15.3）/*ffunction*/

#definex00.0/*initx0,x1*/

#definex11.0

#defineMAXREPT1000/*maxiterationtimes*/

doublex_k=x0,x_k1=x1,x_k2=x1;

Got...%f\n"

x_k2）;

x_k2=x_k1-（f（x_k1）*（x_k1-x_k））/（f（x_k1）-f（x_k））;

if（x_k2-x_k1<

epsilon&

x_k2-x_k1>

-epsilon）

!

Root:

%f\n"

x_k=x_k1;

x_k1=x_k2;

After%drepeat,nosolved.\n"

MAXREPT）;

_

5.运行结果

数值计算方法与算法专业:

5实验项目：

Courant求解方程组指导教师：

一．实验目的：

用C语言实现直接分解法的算法，掌握直接分解法与高斯消去法的不同之处

掌握怎样将已给的算法在计算机上实现，分析算法的优缺点，找到实现直接分解算法最佳的结构；

体会从键盘读入二维数组的源代码，了解计算机内部存储二维数组得规则；

体会Courant分解先分解U的行，再分解L列的算法实现规则；

将下面的方程组用直接分解法实现，并比较其与高斯消元法的不同。

1、用Corant直接分解法求解下列方程组：

aa...axb

aa...ax=b

.....

1.courant的算法描述：

Step1输入矩阵A及列向量b；

Step2将矩阵A分解为A=LU

...

U=...

L=.....

.....

1

Step3对k=1,2,...,n,

=-i=k,k+2,...,n;

=-j=k+1,k+2,...,n

记=y,解=:

=,i=1,2,..,n

Step4再由=得

，j=n,n-1,...,2,1

2.用C语言实现算法的程序代码

intn,i,k;

intj;

intmi;

doublemx,tmp;

staticdoublea[MAX_N][MAX_N],b[MAX_N],x[MAX_N],y[MAX_N];

staticdoublel[MAX_N][MAX_N],u[MAX_N][MAX_N];

\nInputnvalue（dimofAx=b:

TheinputnislargerthanMAX_N,pleaseredefinetheMAX_N.\n"

Nowinputanumbermatrixa（i,j）,i,j=0,...%d:

n-1）;

for（j=0;

j<

j++）

{scanf（"

a[i][j]）;

Nowinputthemaxtrixb（i）,i=0,...,%d:

b[i]）;

u[i][j]=1;

for（k=0;

k<

k++）

for（i=k;

l[i][k]=a[i][k];

=k-1;

l[i][k]-=（l[i][j]*u[j][k]）;

for（i=k+1;

u[k][j]=a[k][j];

u[k][j]-=（l[k][i]*u[i][j]）;

u[k][j]/=l[k][k];

y[i]=b[i];

=i-1;

y[i]-=（l[i][j]*y[j]）;

y[i]/=l[i][i];

for（i=n-1;

i>

=0;

i--）

x[i]=y[i];

for（j=i+1;

x[i]-=（u[i][j]*x[j]）;

Solve...x_i=\n"

x[i]）;

6实验项目：

超松弛迭代求解方程组指导教师：

姓名宋元台学号：

1．实验目的及要求

熟练掌握使用超松弛迭代求解线性方程组

写出超松弛迭代的算法描述并使用C语言实现超松弛代法

2．实验内容

用超松弛迭代（ω作为参数）求解方程组

1.算法描述：

Step1:

输入矩阵A及列向量C；

Step2：

按因子为ω的超松弛迭代公式

......

Step3求解AX=C

#defineMAXREPT100

#defineepslion0.00001

inti,j,k;

doubleerr,w;

staticdoublea[MAX_N][MAX_N],b[MAX_N][MAX_N],c[MAX_N],g[MAX_N];

staticdoublex[MAX_N],nx[MAX_N];

\nInputnvalue（dimofAX=C）:

nowinputthematerixa（i,j）i,j=0...%d\n"

nowinputthemaxtrixc（i）,i=0...%d:

c[i]）;

nowinputthewvalue:

w）;

if（w<

1||w>

=2）

wmustbetween1and2"

b[i][j]=-a[i][j]/a[i][i];

g[i]=c[i]/a[i][i];

nx[j]=g[j];

j;

nx[j]+=b[j][k]*nx[k];

for（k=j+1;

nx[j]+=b[j][k]*x[k];

nx[j]=（1-w）*x[j]+w*nx[j];

err=0;

if（err<

fabs（nx[j]-x[j]））err=fabs（nx[j]-x[j]）;

x[j]=nx[j];

epslion）

Solve...x_i=\n"

After%drepeat,noresult...\n"

re