数值分析实验.docx

1、数值分析实验数值分析实验第二章实验内容：拉格朗日差值问题描述： 利用拉格朗日差值来计算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 #define A 3.14159265

2、4int main() double a11; cout原始的Sin函数的值是：endl; for(int i=0;i11;i+) ai=sin(A*i/20); coutait; coutendl固定节点处的一次拉格朗日差值的值是:endl; double b8; for(int j=0;j8;j+) bj=(aj+aj+1)/2; coutbjt; coutendl固定节点处的二次拉格朗日差值的值是:endl; double c8; for(int s=0;s8;s+) cs=(3*as+6*as+1-as+2)/8; coutcst; coutendl插值节点出的Sin真值是:endl;

3、 double d8; for(int t=0;t8;t+) dt=sin(A*(t+0.5)/20); coutdtt; double f8,h8; coutendl一次拉格朗日插值与真值的误差是:endl; for(int p=0;p8;p+) fp=bp-cp; coutfpt; coutendl二次拉格朗日插值与真值的误差是:endl; for(int q=0;q8;q+) hq=bq-dq; couthqt; coutendl; return 0;测试结果：心得： 在使用拉格朗日差值时，如果区间比较大时，而且不适合用高阶拉格朗日差值时，我们可以先利用样条的思想，把大区间分成m小段，然

4、后在每个小区间上使用低阶的拉格朗日差值，这样得到的结果较之前就精确了很多。从程序运行的结果可以看出这样的结论。第三章实验内容：自动选取步长的龙贝格积分问题描述： 利用龙贝格积分程序，自动选取步长计算下面的三个积分：（1）1.110720735（2）0.500000000（3）算法原码：#include #include #include double F1(double x);/被积分函数void LBG(double a,double b,double real);/龙贝格积分函数可输入积分上下限和函/数真值#define n 100#define PI 3.141592654void ma

5、in() LBG(0,1,1.110720735);/可以修改函数F1的return值来球不同的函数积分/龙呗格积分函数体void LBG(double a,double b,double real) int i,j; double h=b-a; double r3n=0.0; r11=0.5*h*(F1(a)+F1(b); printf(%.7f ,r11);coutendl; for(i=2;i=n;i+) r21=0.5*r11; for(int k=1;k=pow(2,i-2);k+) r21+=F1(a+(k-0.5)*h)*0.5*h; for(j=2;j=i;j+) r2j=r2

6、j-1+(r2j-1-r1j-1)/(pow(4,j-1)-1); printf(%.7f ,r2i); for(j=1;j=i;j+) r1j=r2j; if(fabs(r2i-real)0.000001) break;/满足精度要求时，自动跳出 coutendl; h=h/2; double F1(double x) return (1+x*x)/(1+x*x*x*x);/可修改的函数测试结果：(1) 第一个积分求解情况，真值为1.110720735，控制误差为0.000001.(2)第二个积分求解情况，真值为0.500000000，控制误差为0.000001.(3) 第三个积分求解情况，

7、真值为，控制误差为0.000001.心得： 在求解积分时，由于电脑的字长限制只能表示一部分数据，所以在求解第二个积分的值时，原本精确的解是0.5，但是程序运行的结果是0.4999993.这样的结果虽然不是完全的符合，但是这是电脑可以做的最好的答案了。第四章实验内容：常微分方程解的四阶龙格库塔方法和阿达姆斯方法问题描述： 完成以下功能：(1) 用C+语言描述：阿达姆斯预报校正式解常微分方程初值问题的过程；(2) 调用龙格库塔过程求解下面的一介常微分方程前四点的函数值，并用这些值作阿达姆斯方法的初值； 0x1 h=0.05(3) 调用阿达姆斯过程，求解常微分方程初值问题；(4) 计算解析解在节点的

8、值；解析解的表达式是： (5) 计算解析解与数值解的差；(6) 比较龙格库塔方法和阿达姆斯方法解的精度。算法源码：#include #include #include #define n 21void LGKT(double h0,double y0);/四阶龙格库塔方法函数void ADMS(double h0,double y0);/阿达姆斯预报校正公式函数double TF(double x);/计算常微分方程精确初值的函数double F(double x,double y);/f(x,y)的值输出函数void main() LGKT(0.05,0.0); ADMS(0.05,0.0)

9、;/四阶龙格库塔方法void LGKT(double h0,double y0) int i; double k1,k2,k3,k4,yyn,h=h0,x,y; yy0=y0; for(i=0;in;i+) x=i*h;y=yyi;k1=F(x,y); x=i*h+0.5*h;y=yyi+0.5*h*k1;k2=F(x,y); x=i*h+0.5*h;y=yyi+0.5*h*k2;k3=F(x,y); x=i*h+h;y=yyi+h*k3;k4=F(x,y); yyi+1=yyi+h*(k1+2*k2+2*k3+k4)/6; cout项t四阶龙格库塔公式值t常微分方程的精确值t与解析解的误差e

10、ndl; for(i=1;in;i+) printf(yy%d %.10f %.10f %.10f;,i,yyi,TF(h0*i),TF(h0*i)-yyi); coutendl; /阿达姆斯预报校正方法void ADMS(double h0,double y0) int i; double k1,k2,k3,k4,yyn,h=h0,ffn,x,y; yy0=y0; for(i=0;in;i+) if(i3) /用四阶的龙格库塔方法提供初值 x=i*h;y=yyi;k1=F(x,y),ffi=k1; x=i*h+0.5*h;y=yyi+0.5*h*k1;k2=F(x,y); x=i*h+0.5

11、*h;y=yyi+0.5*h*k2;k3=F(x,y); x=i*h+h;y=yyi+h*k3;k4=F(x,y); yyi+1=yyi+h*(k1+2*k2+2*k3+k4)/6; else /用阿达姆斯校正公式求解 x=i*h;y=yyi;ffi=F(x,y); yyi+1=yyi+h*(55*ffi-59*ffi-1+37*ffi-2-9*ffi-3)/24; x=i*h+h;y=yyi+1;ffi+1=F(x,y); yyi+1=yyi+h*(9*ffi+1+19*ffi-5*ffi-1+ffi-2)/24; cout项t阿达姆斯预报校正公式值t常微分方程的精确值t与解析解的误差end

12、l; for(i=1;in;i+) printf(yy%d %.10f %.10f %.10f;,i,yyi,TF(h0*i),TF(h0*i)-yyi); coutendl; double F(double x,double y) return (2*x*x*x-2*x*y);/可以更改的f(x,y)double TF(double x) return (exp(-1)*x*x)+x*x-1);/可以更改的精确常微分方程值运行结果分析：心得：在用计算机求解数学问题是，要找到适合计算机做的方法。不一定要拘谨于固定的格式，只要是可以利用到已经求解的结果得到更精确的解就可以。第五章实验内容：最小二

13、乘逼近问题描述： 按照下表所示的数据求实验数据的线性最小二乘多项式： XiYiXiYiXiYiXiYi0.04026.50.04128.10.05525.20.05626.00.06224.00.07125.00.07426.40.07426.40.07827.20.08225.60.09025.00.09226.80.10024.80.10527.00.12025.00.12326.9XiYiXiYi0.13026.90.14026.2算法源码：#include #include #define N 17void ZXEC(double x,double y);/最小二乘法的函数void m

14、ain() double xN,yN;/二个数组用于存放最小二乘的x和y值 x0=0.040;y0=26.5;x1=0.041;y1=28.1; x2=0.055;y2=25.2;x3=0.056;y3=26.0; x4=0.062;y4=24.0;x5=0.071;y5=25.0; x6=0.074;y6=26.4;x7=0.078;y7=27.2; x8=0.082;y8=25.6;x9=0.090;y9=25.0; x10=0.092;y10=26.8;x11=0.100;y11=24.8; x12=0.105;y12=27.0;x13=0.120;y13=25.0; x14=0.123

15、;y14=26.9;x15=0.130;y15=26.9; x16=0.140;y16=26.2; ZXEC(x,y);/最小二乘法的函数体void ZXEC(double xN,double yN) double x1=0.0,x2=0.0,x8=0.0,y1=0.0,y8=0.0,xy=0.0,a=0.0,b=0.0; for(int i=0;iN;i+)/计算Xi,Yi,XiYi三项的求和 x1+=xi;x2+=xi*xi; y1+=yi; xy+=xi*yi; x8=x1/N;y8=y1/N;/x和y的均值 b=(xy-N*x8*y8)/(x2-N*x8*x8);/线性表达是的斜率b

16、a=y8-b*x8;/截距a couty=a+(b)xendl;运行结果：心得： 在求解低阶的最小二乘多项式时，可以将低阶的方程组化简后再去求解方程的系数。在本题中就是将未知量化简后再用给出的数据求解。第六章实验内容：非线性方程的求根问题描述：用程序解决下列问题：(1) 用二分法写一个求函数飞f(x)=0实根的标准过程，并解非线性方程在x=0附近的根。根的近似值是0.3604217(2) 用弦截法写一个标准过程并解非线性方程，在4.0附近的根。根的近似值是3.73307903.(3) 用二分法或弦截法写一个标准过程，求多项式.在1.0附近的根。根的近似值是1.07816259.(4) 设计一个

17、不用除法操作，计算的牛顿迭代格式并写出标准过程，其中。取计算并输出迭代次数，。算法源码：#include #include double F(double x);void EFF(double ,double ,double ,double );/二分法求解函数void XJF(double ,double ,double ,double );/弦截法求解函数void Newton(double ,double );/牛顿法求倒数#define M 1000000void main() EFF(0.0,1.0,0.000005,0.3604217); XJF(1.0,1.5,0.000005,

18、1.07816259); cout除法a=0.1,取初值为h=8.8。endl; Newton(8.8,0.1); cout除法a=0.2,取初值为h=4.0。endl; Newton(4.0,0.2); cout除法a=0.5,取初值为h=3.0。endl; Newton(3.0,0.5); cout除法a=0.8,取初值为h=1.0。endl; Newton(1.0,0.8); cout除法a=2.5,取初值为h=0.3。endl; Newton(0.3,2.5); cout除法a=10,取初值为h=0.15。endl; Newton(0.15,10);/二分法求解过程的函数体void E

19、FF(double a,double b,double eps,double real) double x0=a,x1=b, result=0.0,x=(b-a)/2; for(int i=1;iM;i+) result=F(x); if(real-eps)x&x(real+eps) cout二分法求解的方程的根为：xendl; break; if(result0) x0=x; else x1=x; x=(x1-x0)/2+x0; /弦截法求解过程的函数体void XJF(double a,double b,double eps,double real) double x0=a,x1=b,y0

20、=F(a),y1=F(b); for(int i=1;iM;i+) double temp0=x1,temp1=y1; x1=x1-y1*(x1-x0)/(y1-y0); if(real-eps)x1&x1(real+eps) cout弦截法求解的方程的根为：x1endl; break; y1=F(x1); x0=temp0;y0=temp1; /牛顿法求解倒数的函数体void Newton(double h0,double a) double h100; h0=h0; for(int i=0;iM;i+) hi+1=hi*(2-a*hi); couti+1thi+1endl; if(hi+1

21、-hi)0) double e=(-1)*(hi+1-hi); if(e0.001)break; else if(hi+1-hi)0.001) break; /在每个小题中这个函数的内部是不同的,在这里是不能在同一个程序中求解所有问题的，要逐个来求解double F(double x) 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) 用二分法或弦截法写一

22、个标准过程，求多项式.在1.0附近的根。根的近似值是1.07816259.(4) 设计一个不用除法操作，计算的牛顿迭代格式并写出标准过程，其中。取计算并输出迭代次数，。心得： 在解决计算机上解决实际问题时，例如求除法，我们可以根据一些数学的算法，来转化一些计算机很难解决的问题。 在利用二分法或弦截法解决非线性方程求根问题时，我们发现选取很好的初值点式很重要的。如果初值选取的不好是，我们的收敛与否和收敛速度都有一定的影响。第七章实验内容：追赶法求解线性方程组问题问题描述： 用追赶法解三对角线方程组。(1) 写出用追赶法解三对角线方程组的标准过程。(2) 建立方程组的系数矩阵和自由项 置未知量个数

23、 置 置。(3) 将未知向量置零，调用追赶法标准过程求解。(4) 分析计算结果精度。算法源码：#include #include #define M 40void ZGF(double aMM,double bM);/追赶法求解函数void main() /将方程组的系数矩阵aMM和bM,求解出来！ int i,j; double aMM=0.0,fM=0.0; for(i=0;iM;i+) aii=2.0; for(i=0;iM-1;i+) ai+1i=aii+1=-1.0; for(i=1;iM+1;i+) for(j=1;jM+1;j+) fi-1=fi-1+2*j*ai-1j-1; Z

24、GF(a,f);/追赶法的函数体void ZGF(double aMM,double bM) double lMM=0.0,uMM=0.0,xM=0.0; int i; l00=a00; u01=a01/l00; for( i=1;iM;i+) lii-1=aii-1; lii=aii-lii-1*ui-1i; uii+1=aii+1/lii; double yM; y0=b0/l00; for( i=1;i=0;i-) xi=yi-uii+1*xi+1; for(i=0;iM;i+) coutxit;运行结果：心得： 追赶法是一种求解特殊系数矩阵（三对角线）方程组的方法，从中我们可以发现，这样的解法比原始的解法要节约存储空间和运算量。这提醒我们在用计算机解数学上的一些问题时，对一些特殊的情况可以寻找更简便的解法。