数值分析实验.docx

数值分析实验.docx

《数值分析实验.docx》由会员分享，可在线阅读，更多相关《数值分析实验.docx（22页珍藏版）》请在冰豆网上搜索。

数值分析实验

数值分析实验

第二章

实验内容：

拉格朗日差值

问题描述：

利用拉格朗日差值来计算Sin函数的值，并比较插值的值与真实值的误差。

（1）计算函数表sinXi,其中Xi=i*3.1415/20,i=0,1,……,10。

（2）用一次和二次拉格朗日差值公式对函数表进行差值得L1（Xj）,L2（Xj），其中Xj=3.1415*（j+0.5）/20,j=0,1,…….,8。

（3）计算函数值Yi=sin（Xj）,j=0,1,……,8。

（4）计算Ei

（1）=yi-L1（Xj）;Ei

（2）=yi-L2（Xj）;i=0,1,….,8.

原码算法：

#include

#include

#defineA3.141592654

intmain（）

{

doublea[11];

cout<<"原始的Sin函数的值是：

"<

for（inti=0;i<11;i++）

{

a[i]=sin（A*i/20）;

cout<

}

cout<

"<

doubleb[8];

for（intj=0;j<8;j++）

{

b[j]=（a[j]+a[j+1]）/2;

cout<

}

cout<

"<

doublec[8];

for（ints=0;s<8;s++）

{

c[s]=（3*a[s]+6*a[s+1]-a[s+2]）/8;

cout<

}

cout<

"<

doubled[8];

for（intt=0;t<8;t++）

{

d[t]=sin（A*（t+0.5）/20）;

cout<

}

doublef[8],h[8];

cout<

"<

for（intp=0;p<8;p++）

{

f[p]=b[p]-c[p];

cout<

}

cout<

"<

for（intq=0;q<8;q++）

{

h[q]=b[q]-d[q];

cout<

}

cout<

return0;

}

测试结果：

心得：

在使用拉格朗日差值时，如果区间比较大时，而且不适合用高阶拉格朗日差值时，我们可以先利用样条的思想，把大区间分成m小段，然后在每个小区间上使用低阶的拉格朗日差值，这样得到的结果较之前就精确了很多。

从程序运行的结果可以看出这样的结论。

第三章

实验内容：

自动选取步长的龙贝格积分

问题描述：

利用龙贝格积分程序，自动选取步长计算下面的三个积分：

（1）

1.110720735

（2）

0.500000000

（3）

算法原码：

#include

#include

#include

doubleF1（doublex）;//被积分函数

voidLBG（doublea,doubleb,doublereal）;//龙贝格积分函数可输入积分上下限和函//数真值

#definen100

#definePI3.141592654

voidmain（）

{

LBG（0,1,1.110720735）;//可以修改函数F1的return值来球不同的函数积分

}

//龙呗格积分函数体

voidLBG（doublea,doubleb,doublereal）

{

inti,j;

doubleh=b-a;

doubler[3][n]={0.0};

r[1][1]=0.5*h*（F1（a）+F1（b））;

printf（"%.7f",r[1][1]）;cout<

for（i=2;i<=n;i++）

{

r[2][1]=0.5*r[1][1];

for（intk=1;k<=pow（2,i-2）;k++）

r[2][1]+=F1（a+（k-0.5）*h）*0.5*h;

for（j=2;j<=i;j++）

r[2][j]=r[2][j-1]+（r[2][j-1]-r[1][j-1]）/（pow（4,j-1）-1）;

printf（"%.7f",r[2][i]）;

for（j=1;j<=i;j++）

{

r[1][j]=r[2][j];

}

if（fabs（r[2][i]-real）<0.000001）break;//满足精度要求时，自动跳出

cout<

h=h/2;

}

}

doubleF1（doublex）

{

return（1+x*x）/（1+x*x*x*x）;//可修改的函数

}

测试结果：

（1）第一个积分求解情况，真值为1.110720735，控制误差为0.000001.

（2）第二个积分求解情况，真值为0.500000000，控制误差为0.000001.

（3）第三个积分求解情况，真值为

，控制误差为0.000001.

心得：

在求解积分时，由于电脑的字长限制只能表示一部分数据，所以在求解第二个积分的值时，原本精确的解是0.5，但是程序运行的结果是0.4999993.这样的结果虽然不是完全的符合，但是这是电脑可以做的最好的答案了。

第四章

实验内容：

常微分方程解的四阶龙格库塔方法和阿达姆斯方法

问题描述：

完成以下功能：

（1）用C++语言描述：

阿达姆斯预报校正式解常微分方程初值问题的过程；

（2）调用龙格库塔过程求解下面的一介常微分方程前四点的函数值，并用这些值作阿达姆斯方法的初值；

0

（3）调用阿达姆斯过程，求解常微分方程初值问题；

（4）计算解析解在节点的值；

解析解的表达式是：

（5）计算解析解与数值解的差；

（6）比较龙格库塔方法和阿达姆斯方法解的精度。

算法源码：

#include

#include

#include

#definen21

voidLGKT（doubleh0,doubley0）;//四阶龙格——库塔方法函数

voidADMS（doubleh0,doubley0）;//阿达姆斯预报校正公式函数

doubleTF（doublex）;//计算常微分方程精确初值的函数

doubleF（doublex,doubley）;//f（x,y）的值输出函数

voidmain（）

{

LGKT（0.05,0.0）;

ADMS（0.05,0.0）;

}

//四阶龙格——库塔方法

voidLGKT（doubleh0,doubley0）

{

inti;

doublek1,k2,k3,k4,yy[n],h=h0,x,y;

yy[0]=y0;

for（i=0;i

{

x=i*h;y=yy[i];k1=F（x,y）;

x=i*h+0.5*h;y=yy[i]+0.5*h*k1;k2=F（x,y）;

x=i*h+0.5*h;y=yy[i]+0.5*h*k2;k3=F（x,y）;

x=i*h+h;y=yy[i]+h*k3;k4=F（x,y）;

yy[i+1]=yy[i]+h*（k1+2*k2+2*k3+k4）/6;

}

cout<<"项"<<"\t"<<"四阶龙格库塔公式值"<<"\t"<<"常微分方程的精确值"<<"\t"<<"与解析解的误差"<

for（i=1;i

{

printf（"yy[%d]%.10f%.10f%.10f;",i,yy[i],TF（h0*i）,TF（h0*i）-yy[i]）;

cout<

}

}

//阿达姆斯预报——校正方法

voidADMS（doubleh0,doubley0）

{

inti;

doublek1,k2,k3,k4,yy[n],h=h0,ff[n],x,y;

yy[0]=y0;

for（i=0;i

{

if（i<3）

{

//用四阶的龙格库塔方法提供初值

x=i*h;y=yy[i];k1=F（x,y）,ff[i]=k1;

x=i*h+0.5*h;y=yy[i]+0.5*h*k1;k2=F（x,y）;

x=i*h+0.5*h;y=yy[i]+0.5*h*k2;k3=F（x,y）;

x=i*h+h;y=yy[i]+h*k3;k4=F（x,y）;

yy[i+1]=yy[i]+h*（k1+2*k2+2*k3+k4）/6;

}

else

{

//用阿达姆斯校正公式求解

x=i*h;y=yy[i];ff[i]=F（x,y）;

yy[i+1]=yy[i]+h*（55*ff[i]-59*ff[i-1]+37*ff[i-2]-9*ff[i-3]）/24;

x=i*h+h;y=yy[i+1];ff[i+1]=F（x,y）;

yy[i+1]=yy[i]+h*（9*ff[i+1]+19*ff[i]-5*ff[i-1]+ff[i-2]）/24;

}

}

cout<<"项"<<"\t"<<"阿达姆斯预报校正公式值"<<"\t"<<"常微分方程的精确值"<<"\t"<<"与解析解的误差"<

for（i=1;i

{

printf（"yy[%d]%.10f%.10f%.10f;",i,yy[i],TF（h0*i）,TF（h0*i）-yy[i]）;

cout<

}

}

doubleF（doublex,doubley）

{

return（2*x*x*x-2*x*y）;//可以更改的f（x,y）

}

doubleTF（doublex）

{

return（exp（（-1）*x*x）+x*x-1）;//可以更改的精确常微分方程值

}

运行结果分析：

心得：

在用计算机求解数学问题是，要找到适合计算机做的方法。

不一定要拘谨于固定的格式，只要是可以利用到已经求解的结果得到更精确的解就可以。

第五章

实验内容：

最小二乘逼近

问题描述：

按照下表所示的数据求实验数据的线性最小二乘多项式：

Xi

Yi

Xi

Yi

Xi

Yi

Xi

Yi

0.040

26.5

0.041

28.1

0.055

25.2

0.056

26.0

0.062

24.0

0.071

25.0

0.074

26.4

0.074

26.4

0.078

27.2

0.082

25.6

0.090

25.0

0.092

26.8

0.100

24.8

0.105

27.0

0.120

25.0

0.123

26.9

Xi

Yi

Xi

Yi

0.130

26.9

0.140

26.2

算法源码：

#include

#include

#defineN17

voidZXEC（doublex[],doubley[]）;//最小二乘法的函数

voidmain（）

{

doublex[N],y[N];//二个数组用于存放最小二乘的x和y值

x[0]=0.040;y[0]=26.5;x[1]=0.041;y[1]=28.1;

x[2]=0.055;y[2]=25.2;x[3]=0.056;y[3]=26.0;

x[4]=0.062;y[4]=24.0;x[5]=0.071;y[5]=25.0;

x[6]=0.074;y[6]=26.4;x[7]=0.078;y[7]=27.2;

x[8]=0.082;y[8]=25.6;x[9]=0.090;y[9]=25.0;

x[10]=0.092;y[10]=26.8;x[11]=0.100;y[11]=24.8;

x[12]=0.105;y[12]=27.0;x[13]=0.120;y[13]=25.0;

x[14]=0.123;y[14]=26.9;x[15]=0.130;y[15]=26.9;

x[16]=0.140;y[16]=26.2;

ZXEC（x,y）;

}

//最小二乘法的函数体

voidZXEC（doublex[N],doubley[N]）

{

doublex1=0.0,x2=0.0,x8=0.0,y1=0.0,y8=0.0,xy=0.0,a=0.0,b=0.0;

for（inti=0;i

{

x1+=x[i];x2+=x[i]*x[i];

y1+=y[i];

xy+=x[i]*y[i];

}

x8=x1/N;y8=y1/N;//x和y的均值

b=（xy-N*x8*y8）/（x2-N*x8*x8）;//线性表达是的斜率b

a=y8-b*x8;//截距a

cout<<"y="<

}

运行结果：

心得：

在求解低阶的最小二乘多项式时，可以将低阶的方程组化简后再去求解方程的系数。

在本题中就是将未知量化简后再用给出的数据求解。

第六章

实验内容：

非线性方程的求根

问题描述：

用程序解决下列问题：

（1）用二分法写一个求函数飞f（x）=0实根的标准过程，并解非线性方程

在x=0附近的根。

根的近似值是0.3604217

（2）用弦截法写一个标准过程并解非线性方程

，在4.0附近的根。

根的近似值是3.73307903.

（3）用二分法或弦截法写一个标准过程，求多项式

.在1.0附近的根。

根的近似值是1.07816259.

（4）设计一个不用除法操作，计算

的牛顿迭代格式并写出标准过程，其中

。

取

计算并输出迭代次数，

。

算法源码：

#include

#include

doubleF（doublex）;

voidEFF（double,double,double,double）;//二分法求解函数

voidXJF（double,double,double,double）;//弦截法求解函数

voidNewton（double,double）;//牛顿法求倒数

#defineM1000000

voidmain（）

{

EFF（0.0,1.0,0.000005,0.3604217）;

XJF（1.0,1.5,0.000005,1.07816259）;

cout<<"除法a=0.1,取初值为h=8.8。

"<

Newton（8.8,0.1）;

cout<<"除法a=0.2,取初值为h=4.0。

"<

Newton（4.0,0.2）;

cout<<"除法a=0.5,取初值为h=3.0。

"<

Newton（3.0,0.5）;

cout<<"除法a=0.8,取初值为h=1.0。

"<

Newton（1.0,0.8）;

cout<<"除法a=2.5,取初值为h=0.3。

"<

Newton（0.3,2.5）;

cout<<"除法a=10,取初值为h=0.15。

"<

Newton（0.15,10）;

}

//二分法求解过程的函数体

voidEFF（doublea,doubleb,doubleeps,doublereal）

{

doublex0=a,x1=b,result=0.0,x=（b-a）/2;

for（inti=1;i

{

result=F（x）;

if（（real-eps）

{

cout<<"二分法求解的方程的根为：

"<

break;

}

if（result<0）x0=x;

elsex1=x;

x=（x1-x0）/2+x0;

}

}

//弦截法求解过程的函数体

voidXJF（doublea,doubleb,doubleeps,doublereal）

{

doublex0=a,x1=b,y0=F（a）,y1=F（b）;

for（inti=1;i

{

doubletemp0=x1,temp1=y1;

x1=x1-y1*（x1-x0）/（y1-y0）;

if（（real-eps）

{

cout<<"弦截法求解的方程的根为：

"<

break;

}

y1=F（x1）;

x0=temp0;y0=temp1;

}

}

//牛顿法求解倒数的函数体

voidNewton（doubleh0,doublea）

{

doubleh[100];

h[0]=h0;

for（inti=0;i

{

h[i+1]=h[i]*（2-a*h[i]）;

cout<

if（（h[i+1]-h[i]）<0）

{

doublee=（-1）*（h[i+1]-h[i]）;

if（e<0.001）break;

}

elseif（（h[i+1]-h[i]）<0.001）break;

}

}

//在每个小题中这个函数的内部是不同的,在这里是不能在同一个程序中求解所有问题的，要逐个来求解

doubleF（doublex）

{

return（2*x*x*x+4*x*x-2*x-5）;

}

运行结果：

（1）用二分法写一个求函数飞f（x）=0实根的标准过程，并解非线性方程

在x=0附近的根。

根的近似值是0.3604217

（2）用弦截法写一个标准过程并解非线性方程

，在4.0附近的根。

根的近似值是3.73307903.

（3）用二分法或弦截法写一个标准过程，求多项式

.在1.0附近的根。

根的近似值是1.07816259.

（4）设计一个不用除法操作，计算

的牛顿迭代格式并写出标准过程，其中

。

取

计算并输出迭代次数，

。

心得：

在解决计算机上解决实际问题时，例如求除法，我们可以根据一些数学的算法，来转化一些计算机很难解决的问题。

在利用二分法或弦截法解决非线性方程求根问题时，我们发现选取很好的初值点式很重要的。

如果初值选取的不好是，我们的收敛与否和收敛速度都有一定的影响。

第七章

实验内容：

追赶法求解线性方程组问题

问题描述：

用追赶法解三对角线方程组。

（1）写出用追赶法解三对角线方程组的标准过程。

（2）建立方程组

的系数矩阵

和自由项

●置未知量个数

●置

●置

。

（3）将未知向量

置零，调用追赶法标准过程求解。

（4）分析计算结果精度。

算法源码：

#include

#include

#defineM40

voidZGF（doublea[M][M],doubleb[M]）;//追赶法求解函数

voidmain（）

{

//将方程组的系数矩阵a[M][M]和b[M],求解出来！

inti,j;

doublea[M][M]={0.0},f[M]={0.0};

for（i=0;i

a[i][i]=2.0;

for（i=0;i

a[i+1][i]=a[i][i+1]=-1.0;

for（i=1;i

for（j=1;j

f[i-1]=f[i-1]+2*j*a[i-1][j-1];

ZGF（a,f）;

}

//追赶法的函数体

voidZGF（doublea[M][M],doubleb[M]）

{

doublel[M][M]={0.0},u[M][M]={0.0},x[M]={0.0};

inti;

l[0][0]=a[0][0];

u[0][1]=a[0][1]/l[0][0];

for（i=1;i

{

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

l[i][i]=a[i][i]-l[i][i-1]*u[i-1][i];

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

}

doubley[M];

y[0]=b[0]/l[0][0];

for（i=1;i

y[i]=（b[i]-l[i][i-1]*y[i-1]）/l[i][i];

x[M-1]=y[M-1];

for（i=M-2;i>=0;i--）

x[i]=y[i]-u[i][i+1]*x[i+1];

for（i=0;i

cout<

}

运行结果：

心得：

追赶法是一种求解特殊系数矩阵（三对角线）方程组的方法，从中我们可以发现，这样的解法比原始的解法要节约存储空间和运算量。

这提醒我们在用计算机解数学上的一些问题时，对一些特殊的情况可以寻找更简便的解法。