5《数值逼近》课程设计报告.docx

1、5数值逼近课程设计报告数值逼近课程设计报告一、目的意义 (1)进一步熟悉掌握复化梯形和复化抛物线公式(2)学会比较复化梯形公式和复化抛物线公式如何达到所要求的精度(3) 提高编程能力(4)通过数值方法求出很难求得原函数的积分和解析表达是没有明确的给出积分的近似值二、内容要求 积分计算问题：分别用复化梯形和复化Simpsom求积公式计算积分，并比较计算量（精度为10-8）。三、问题解决的方法与算法方法：解决该积分问题时，运用了数值积分近似解法的方法，运用复化梯形和复化Simpsom求积公式进行计算3.1 复化梯形积分3.1.1复化梯形积分公式表达式3.1.2复化梯形积分误差表达式3.2复化抛物线

2、积分3.2.1复化抛物线积分公式表达式3.2.2复化抛物线积分误差表达式3.3算法3.3.1复化梯形积分算法第一步：根据精度计算n的值，输入两端点的值，计算步长h第二步：根据步长计算出各个节点xi的值，i=0,1,2,n第三步：根据xi计算出各个节点对应yi的值，i=0,1,2,n第四步：对各个节点的值进行求和第五步：输出最终的积分的值3.3.2复化抛物线积分算法第一步：根据精度计算n的值，端点a,b的值，计算步长h第二步：根据步长计算出各个节点xi的值，i=0,1,2,n第三步：根据xi计算出各个节点对应yi的值，i=0,1,2,n第四步：对各个节点的值进行求和,分情况，对左右端点先求和，对

3、剩下的端点，奇数的求和后乘以4倍，偶数的求和后乘以2倍，最终将各个值进行加和第五步：对加和的值乘以步长除以3第六步：输出最终的积分的值四、计算程序/ 复化梯形公式.cpp : 定义控制台应用程序的入口点。/n+1点的复化梯形公式#include/#include stdafx.h#include#include # include using namespace std;const int ARRAY_LEN (10000);class Comtprotected:double a,b,h,n,xARRAY_LEN,yARRAY_LEN;char fARRAY_LEN;public:void

4、getab() coutab; coutendl;void cal_nh() int c; coutc; coutendl; n=c-1; coutn的值为：nendl; h=(b-a)/n; couth的值为：hendl;void cal_x() int i=0; double temp=0; for(i;in+1;i+) temp=i*h; xi=a+temp; /*coutx的值为：endl; for(i=0;in+1;i+) coutxi=xi ; coutendlendl;*/void get_f() char tempARRAY_LEN; couttemp; strcpy(f,te

5、mp); coutendl;void cal_y() int i=0; double temp=0; for(i=0;in+1;i+) temp=13*(xi-xi*xi)*exp(-1.5*xi); yi=temp; /*couty的值为：endl; for(i=0;in+1;i+) coutyi=yi ; coutendlendl;*/double result() double temp=y0,sum; int i=1; for(i=1;in+1;i+) if(i=n) temp=temp+yi;/当有判断条件时，要先进行判断，不满足时才进行原始计算 else temp=temp+2*y

6、i; sum=h*temp/2; return sum;void display() double m; m=result(); cout积分区间为a b,被积函数为f的积分endl; cout用n+1点复化梯形公式endl计算结果为:setprecision(8)endlmendl; ;int main() Comt com; com.getab(); com.get_f(); com.cal_nh(); com.cal_x(); com.cal_y(); com.display(); return 0;/复化抛物线公式程序#include #include#include# include

7、 using namespace std;const int ARRAY_LEN (10000);class Comsimprotected:double a,b,h,n,xARRAY_LEN,yARRAY_LEN;char fARRAY_LEN;public:void getab() coutab; coutendl;void cal_nh() int c; coutc; n=2*c-1; coutn的值为：n; h=(b-a)/(2*n); couth的值为：h;void cal_x() int i=0; double temp=0; for(i;i2*n+1;i+) temp=i*h;

8、xi=a+temp; /*coutx的值为：endl; for(i=0;i2*n+1;i+) coutxi=xi ; if(i%5=0) coutendl; coutendlendl;*/void get_f() char tempARRAY_LEN; couttemp; strcpy(f,temp);void cal_y() int i=0; double temp=0; for(i=0;i2*n+1;i+) temp=13*(xi-xi*xi)*exp(-1.5*xi); yi=temp; /*couty的值为：endl; for(i=0;i2*n+1;i+) coutyi=setprec

9、ision(8)yi ; if(i%5=0) coutendl; coutendlendl;*/double result() double temp1=y0,temp2=0,temp3=0,sum1=0,sum=0; int i=1; for(i=1;i2*n+1;i+) if(i=2*n) temp1=temp1+yi;/当有判断条件时，要先进行判断，不满足时才进行原始计算 else if(i%2=0) temp2=temp2+2*yi; else temp3=temp3+4*yi; sum1=temp1+temp2+temp3; sum=h*sum1/3; return sum;void

10、 display() double m; m=result(); cout积分区间为a b,被积函数为f的积分endl; cout用n+1点复化抛物线公式endl计算结果为:setprecision(8)m=0&x=52.36*17&x=1884.96)则用分段线性插值进行计算，如果两种情况均不是则输出输入错误第四步：计算出柱高S第五步：输出S四、计算程序4.1用牛顿插值和分段线性插值的程序/ nt.cpp : 定义控制台应用程序的入口点。/#include StdAfx.h#include#define ARRAY_LEN 18double dq(double a) int i=0; int

11、 j=0; double mul=1; double h=52.36; double xARRAY_LEN=0; double yARRAY_LEN=502.75,520.96,525,523.6,514.3,492,451,394.6,326.5,256.7,188.6,132.1,92.2,68.9,59.6,58.2,62.24; double temp=0; for(i=0;iARRAY_LEN;i+) temp=i*h; xi=temp; printf(请输出所有x的值：n); for(i=0;iARRAY_LEN;i+) printf(%f ,xi); if(i%3=0) prin

12、tf(n); printf(n); printf(请输出所有y的值：n); for(j=0;jARRAY_LEN;j+) printf(%lf ,yj); if(j%3=0) printf(n); printf(n); /求差商 /printf(各阶插商值为：n); for(i=0;ii;j-) yj=(yj-yj-1)/(xj-xj-1-i); /求和 double s=y0; for(j=1;jARRAY_LEN;j+) mul=1; for(i=0;i=j-1;i+) mul=mul*(a-xi); s=s+mul*yj; return s;double s(double n) int

13、i=0; int j=0; double s=0; double x2=890.12,1884.96; double y2=80.45,502.75; printf(请输出所有x的值：n); for(i=0;i2;i+) printf(%f ,xi); if(i%2=0) printf(n); printf(n); printf(请输出所有y的值：n); for(j=0;j2;j+) printf(%f ,yi); if(i%2=0) printf(n); /线形函数 for(i=1;i=0&N=52.36*17&N=1884.96) S=s(N); printf(n); printf(x=%

14、f时的柱高的近似值为：%lfn,N,S); else printf(输入的值不能进行计算！n); return 0;4.2用多项式拟合的matlab程序4.2.1前曲线段拟合五次多项式的值clear all;clc;x=0:52.36:16*52.36;y=502.75 520.96 525 523.6 514.3 492 451 394.6 326.5 256.7 188.6 132.1 92.2 68.9 59.6 58.2 62.24;b=polyfit(x,y,5); % 拟合出的五次函数的系数fprintf(运行结果为：nn)%输出语句disp(拟合出的五次多项式P4(x)为：)%输

15、出语句f=poly2str(b,x)%将拟合后的多项式系数（双精度数组）转换为字符形式的函数poly2sym(b);%将该向量转换为多项式xx=linspace(min(x),max(x); % 绘图用到的点的横坐标yy=polyval(b,xx); % 拟合曲线的纵坐标%subplot(2,2,2);plot(x,y,m.,xx,yy,-c); % 绘图，原始数据+拟合曲线xlabel(圆周上任一点x);ylabel(对应的柱高);legend(原始数据,拟合曲线); % 图示title(五次多项式拟合曲线);hold on;fprintf(x=53时柱高的近似值为：n)e=53 104 2

16、08 409m=polyval(b,e)%用于对已经拟合后的多项式系数， %当给出某个点时求其函数值;计算插值多项式在pi/8处的值plot(e,m,b*)4.2.2线性拟合后半段x=52.36*17 1884.96;y=80.45 502.75;b=polyfit(x,y,1); % 拟合出的1次函数的系数f=poly2str(b,x)%将拟合后的多项式系数（双精度数组）转换为字符形式的函数poly2sym(b);%将该向量转换为多项式xx=linspace(min(x),max(x); % 绘图用到的点的横坐标yy=polyval(b,xx); % 拟合曲线的纵坐标plot(x,y,m.,

17、xx,yy,-c); % 绘图，原始数据+拟合曲线xlabel(圆周上任一点x);ylabel(对应的柱高);legend(原始数据,拟合曲线); % 图示title(一次线性拟合曲线);hold on;x=52.36*17 52.36*21 52.36*29 52.36*35fprintf(x对应柱高的近似值为：n)m=polyval(b,e)%用于对已经拟合后的多项式系数 %当给出某个点时求其函数值plot(e,m,b*)五、计算结果与分析5.1用牛顿插值和分段线性插值的运行结果图5.1.1时程序的运行结果图5.1.2时程序的运行结果5.2用多项式拟合的matlab程序运行结果5.2.1前

18、曲线段拟合五次多项式的结果下面为自己给定的x值经过五次多项式拟合后计算的结果，图5.2.1为拟合的曲线图。运行结果为：拟合出的五次多项式P5(x)为：f =-1.678e-011 x5 + 3.6684e-008 x4 - 2.4645e-005 x3 + 0.0041844 x2 - 0.042005 x + 506.7499x=53 104 208 409时柱高的近似值为：m = 512.8908 524.0042 519.3960 337.8354图5.2.1五次多项式拟合图5.2.2线性拟合后半段结果下面为自己给定的x值经过线性拟合后计算的结果，图5.2.2为拟合的结果图。f = 0.

19、42449 x - 297.3974x = 1.0e+003 * 0.8901 1.0996 1.5184 1.8326x对应柱高的近似值为：m = 80.4500 169.3553 347.1658 480.5237图5.2.2一次线性拟合结果图5.3结果分析在求解该题的过程中，我分别运用了C语言和matlab两个软件进行编程实现，在C语言中，在时，运用了牛顿插值多项式来进行逼近原曲线，为了验证程序的准确性，在这里取x=53来进行计算，计算结果为512.241797在520.96 525之间，且接近520.96，即证明了程序的准确性。同理，可以取任意来验证，在这里取x=1880进行验证，运行结果见图5.1.2。在matlab中，进行了不同区间的多项式拟合，拟合的结果见图5.2.1和图5.2.2，并且拟合的多项式表达式见上面的运行结果。比较C语言和matlab程序可知，两者的结果较为接近，但是在matlab中可以得到拟合的多项式的表达式，并且可以可到结果图，让人更直观的看出结果。六、