数值计算上机作业.docx

西 南 交 通 大 学数值计算课程上机作业组号 第1组 组员 学号 班级 任课教师 2017年12月10日西南交通大学数值计算课程上机作业目 录第1题错误!未定义书签。1.1 (a)错误!未定义书签。1.2 (b)错误!未定义书签。第2题12.1 (a)错误!未定义书签。2.2 (b)错误!未定义书签。第3题1第4题24.1 (a)错误!未定义书签。4.2 (b)错误!未定义书签。4.3 (c)错误!未定义书签。4.4 (d)错误!未定义书签。4.5 (e)错误!未定义书签。第5题错误!未定义书签。5.1 (a)错误!未定义书签。5.2 (b)错误!未定义书签。第6题 数值积分26.1 (a)错误!未定义书签。6.2 (b)错误!未定义书签。6.3 (c)错误!未定义书签。附录-程序清单14第II页西南交通大学机械原理课设计程项目进展报告第1题题中所给非线性方程的曲线如下图1所示。由图中可以看出方程根的个数为5和每个根所在大致区间，之后采用截线法求得其根的近似解。图 1截线法迭代计算结果：a）采用截线法求该方程的5个根的计算结果如下：x1=-1.38129848x2=-0.66666666x3=0.20518292x4=0.50000000x5=1.17611556收敛速度分析b）讨论根的收敛速度当f(x)在根的某个区域内一阶导数不为零，二阶导数连续时，截线法是超线性收敛的。此处fx=324x5+225x4-408x3-207x2+70x+16fx=1620x4+900x3-1224x2-414x+70可以看出fx在x1、x2、x3、x4、x5附近都是连续的，且fx1、fx3、fx4、fx5都不为0，所以截线法在根x1、x3、x4、x5处式超线性收敛的，x2处是线性收敛的。第2题 解题要求：取一根附近初值，用Newton法求f的根，并给出前50步迭代结果。其中，f（x）=(1-34x)13（3-1）。结果分析。2.1 (a)表3-1 Newton法迭代结果迭代次数n结果x迭代次数n结果xn=1x=0.64n=1x=0.75n=2x=-1.#INDn=2x=0.75.n=50x=-1.#INDn=50x=0.752.2 (b)第3题 解题要求：用Jacobi方法求该方程组（n=1000）的解。精确到小数点后三位需要迭代步数。计算解的误差（用无穷大范数）。 表4-1n值迭代步数解的精度解的误差5170.000974419030.9970767410580.000944671670.010956971002600.00099923950.9465379110002540.000996285011第4题结果4.1（a）f =- 0.000735246*t9 + 0.0359507*t8 - 0.750705*t7 + 8.74608*t6 - 62.32*t5 + 279.926*t4 - 786.973*t3 + 1323.62*t2 - 1192.52*t + 497.288Y=f（X）图像 插值多项式图像 4.2(b)估计2010年产量：syms t f %定义符号变量 f =- 0.000735246*t9 + 0.0359507*t8 - 0.750705*t7 + 8.74608*t6 - 62.32*t5 + 279.926*t4 - 786.973*t3 + 1323.62*t2 - 1192.52*t + 497.288 %定义函数表达式 subs(f,t,18) ans =-4.3447e+064.3(c)存在龙格现象。对应该列，插值多项式不是一个好的模型，因为当n较大时，该函数在等距节点下的高次插值多项式会发生激烈的震荡。4.4（d）前四组数据插值多项式为： xi=1,2,3,4;yi=67.052,68.008,69.803,72.024;f=lagranz(xi,yi) f =- 0.0688333*t3 + 0.8325*t2 - 1.05967*t + 67.3481998年产量： syms t f %定义符号变量 f=-0.0688333*t3+0.8325*t2-1.05967*t+67.348 %定义函数表达式 subs(f,t,5) ans =74.2580故此时不存在龙格现象。4.5（e）拟合程序： x=1 2 3 4 5 6 7 8 9 10;y=67.052 68.008 69.803 72.024 73.400 72.063 74.669 74.487 74.065 76.777 nafit(x,y,1)ans = 0.9693 66.9034 nafit(x,y,2)ans = -0.0724 1.7652 65.3116 nafit(x,y,3)ans = 0.0182 -0.3729 3.1512 63.7490故一次拟合曲线为：y=0.9693x+66.9034 抛物线拟合曲线为： 三次拟合曲线为：（ii）估计2010年产量：将18带入，得y1=84.3508 y2=73.6276 Y3=105.7934（iii）所以 我觉得 最小二乘拟合直线最能表示表格中的数据。第5题 数值积分 解题要求：复合梯形公式/Simpson公式、Romberg公式计算下列积分，达到精度1210-8。给出所需的小区间数量。6.1 (a)精确解：f*=01x3x2+1dx=1201x2x2+1dx2=1201（1-1x2+1）dx2=1-ln22 （6-1） 思路： 从节点数n=1开始使用对应的求积公式，并每次增加一个节点，当f-f*eps,达到目标精度。6.1.1 复合梯形公式结果 表6-1 复合梯形公式结果表节点个数n步长h结果fn=4056h=0.00024654832f=0.153426414786n=4057h=0.00024648755f=0.153426414783.n=4081h=0.00024503798f=0.153426414724n=4082h=0.00024497795f=0.153426414721n=4083h=0.00024491795f=0.153426414719结果：f=0.153426414719复合梯形公式所需小区间为n=4083。6.1.2 Simpson公式结果表6-2 Simpson公式结果表节点个数n步长h结果fn=1h=0.5f=0.15n=2h=0.25f=0.153235294118.n=25h=0.02f=0.153426403048n=26h=0.019230769f=0.153426404017n=27h=0.018518519f=0.153426404817结果：f=0.153426404817Simpson公式所需小区间数量为2n=54。6.1.3 Romberg公式结果表6-3 Romberg公式结果表0.1534460000000.1535070.153426000000.1537520.15342601547310.1534250.1534260.1534260000.1586760.1534150.1534260.1534260.153426000.1750.1532350.1534270.1534260.1534260.15342600.250.150.1534500.1534270.1534260.1534260果：f=0omberg公式所需小区间数量为m=7。6.2 (b)6.2.1 复合梯形公式结果精确解 思路与与（a）相同，但是，由于0ecosxdx（6-2）无法用初等函数计算，因此根据|fn-fn-1|eps,来判断是否达到目标精度。 表6-4 复合梯形公式结果表节点个数n步长h结果fn=1h=3.1415927f=4.84773078623n=2h=1.5707963f=3.99466171991.n=4h=0.78539816f=3.97746388635n=5h=0.62831853f=3.97746326224n=6h=0.52359878f=3.97746326051结果：f=3.97746326051复合梯形公式所需小区间数量为n=6。6.2.2 Simpson公式结果表6-5 Simpson公式结果表节点个数n步长h结果fn=1h=1.5707963f=3.71030536447n=2h=0.78539816f=3.97173127516.n=4h=0.39269908f=3.97746305189n=5h=0.31415927f=3.97746325993n=6h=0.26179939f=3.97746326051结果：f=3.97746326051Simpson公式所需小区间数量2n=12。6.2.3 Romberg公式结果表6-6 Romberg公式结果表3.97746300000003.9774633.9774630000003.9774633.9774633.977463000003.9774633.9774633.9774633.97746300003.9774633.9774633.9774633.9774633.9774630003.9774633.9774633.9774633.9774633.9774633.977463003.9946613.9717313.9778453.9774573.9774633.9774633.97746304.8477303.7103053.9891593.9776653.9774563.9774633.9774633.97746326051结果：f=3.97746326051Romberg公式所需小区间数量为m=8。6.3 (c)精确解 f*=0231x2+4dx=ln(3+2) （6-3） 思路 从节点数n=1开始使用对应的求积公式，并每次增加一个节点，当f-f*eps,达到目标精度。6.3.1 复合梯形公式结果表6-7 复合梯形公式结果表节点个数n步长h结果fn=3263h=0.0010616309f=1.31695789184n=3264h=0.0010613056f=1.31695789184.n=3289h=0.0010532386f=1.31695789192n=3290h=0.0010529184f=1.31695789192n=3291h=0.0010525985f=1.31695789193结果：f=1.31695789193复合梯形公式所需小区间数量为n=3291。6.3.2 Simpson公式结果表6-8 Simpson公式结果表节点个数n步长h结果fn=1h=1.7320508f=1.30588426284n=2h=0.8660254f=1.31671754654.n=15h=0.11547005f=1.3169578894n=16h=0.10825318f=1.31695789111n=17h=0.10188534f=1.31695789236结果：f=1.31695789236Simpson公式所需小区间数量2n=34。6.3.3 Romberg公式结果表6-9 Romberg公式结果表1.3169440000001.3169051.316957000001.3167461.3169571.31695700001.3161121.3169571.3169571.3169570001.3135811.3169561.3169571.3169571.316957001.3041721.3167171.3169721.3169571.3169571.31695701.2990381.3058841.3174391.3169641.3169571.3169571.31695789692结果：f=1.31695789692Romberg公式所需小区间数量为m=7。第7题 解题要求：给定初值条件y(0)=1，求解方程的准确值。在区间0,1上分别取步长0.1,0.05,0.025用Euler法求解数值解，并作图与准确解比较。在区间0,1上分别取步长0.1,0.05,0.025用改进的Euler法求解数值解，并作图与准确解比较。7.1 (a)7.1.1 方程准确解结果7.1.2 Euler法结果表7-1 步长结果表 h=0.1h=0.05h=0.025XYXYXYX=0.0Y=1.00X=0.00Y=1.0000X=0.000Y=1.000000X=0.1Y=1.00X=0.05Y=1.0000X=0.025Y=1.000000.X=0.8Y=1.25X=0.90Y=1.3825X=0.950Y=1.439375X=0.9Y=1.36X=0.95Y=1.4275X=0.975Y=1.463125X=1.0Y=1.45X=1.00Y=1.4750X=1.000Y=1.487500Euler法图形与准确解图形比较：如图7-1 图7-1 图7-2 7.1.3 改进Euler法结果 表7-2 步长结果表 h=0.1h=0.05h=0.025XYXYXYx=0.0y=1.000x=0.00y=1.00000x=0.000y=1.0000000x=0.1y= 1.005x=0.05y= 1.00125x=0.025y=1.0003125.X=0.8Y=1.320X=0.90Y=1.40500X=0.950Y= 1.4512500X=0.9Y=1.405X=0.95Y= 1.45125X=0.975Y= 1.4753125X=1.0Y=1.500X=1.00Y= 1.50000X=1.000Y= 1.5000000改进Euler法图形与准确解图形比较：如图7-27.2 (b)7.2.1 方程准确解结果7.2.2 Euler法结果表7-3 步长结果表 h=0.1h=0.05h=0.025XYXYXYX=0.0Y=1.0000000X=0.00Y=1.0000000X=0.000Y=1.0000000X=0.1Y=1.0000000X=0.05Y=1.0000000X=0.025Y=1.0000000.X=0.8Y=1.1476561X=0.90Y=1.2468407X=0.950Y= 1.3136620X=0.9Y=1.2211061X=0.95Y=1.2973377X=0.975Y= 1.3433015X=1.0Y=1.3200000X=1.00Y=1.3558801X=1.000Y=1 .3752259Euler法图形与准确解图形比较：如图7-3 图7-3 图7-47.2.3 改进Euler法结果 表7-4 步长结果表 h=0.1h=0.05h=0.025XYXYXYX=0.0Y=1.000X=0.00Y=1.00000X=0.000Y=1.0000000X=0.1Y= 1.005X=0.05Y= 1.00125X=0.025Y=1.0003125.X=0.8Y=1.320X=0.90Y=1.40500X=0.950Y= 1.4512500X=0.9Y=1.405X=0.95Y= 1.45125X=0.975Y= 1.4753125X=1.0Y=1.500X=1.00Y= 1.50000X=1.000Y= 1.5000000改进Euler法图形与准确解图形比较：如图7-47.3 (c)7.3.1 方程准确解结果7.3.2 Euler法结果 表7-5 步长结果表 h=0.1h=0.05h=0.025XYXYXYX=0.0Y=1.0000000X=0.00Y=1.0000000X=0.000Y=1.0000000X=0.1Y=1.2000000X=0.05Y=1.1000000X=0.025Y=1.0500000.X=0.8Y=6.7323294X=0.90Y=1 0.949539X=0.950Y= 14.566841X=0.9Y=9.1559680X=0.95Y= 13.029951X=0.975Y= 15.987108X=1.0Y=12.635236X=1.00Y= 15.570792X=1.000Y= 17.565835Euler法图形与准确解图形比较：如图7-5 图7-5 7.3.3 改进Euler法结果 表7-6 步长结果表 h=0.1h=0.05h=0.025XYXYXYx=0.0y=1.0000000x=0.00y=1.0000000x=0.000y=1.0000000x=0.1y= 1.0005000x=0.05y= 1.0000625x=0.025y= 1.0000078.X=0.8Y=1.1875241X=0.90Y= 1.2754926X=0.950Y= 1.3309303X=0.9Y=1.2669770X=0.95Y= 1.3312651X=0.975Y= 1.3621169X=1.0Y=1.3974094X=1.00Y= 1.3960853X=1.000Y= 1.3957338改进Euler法图形与准确解图形比较 附录第2页程序清单1.Matlab出图源代码：clc;clear;x=-1.4:0.001:1.2;y1=-4+16*x+35*x.2-69*x.3-102*x.4+45*x.5+54*x.6;plot(x,y1);grid on;C+计算源代码：#include#include#includeusing namespace std;double x1,x2,x,e;int i;double y(double s);void main()/求根算法，依次求解5个根for(i=0;ipow(10.0,-15)x=x2-(y(x2)/(y(x2)-y(x1)*(x2-x1);/迭代公式x1=x2;x2=x;e=fabs(x2-x1)/x2);/精度计算公式/输出格式控制，保留到小数点后8位coutfixed;cout.precision(8);coutxi+1=x2endl;/题中所给方程double y(double s)double a;a=54.0*pow(s,6)+45.0*pow(s,5)-102.0*pow(s,4)-69.0*pow(s,3)+35.0*pow(s,2)+16.0*s-4.0;return a;程序清单2.#include#include#include /控制输出有效数字位数函数所在的库using namespace std;double fname(double x) double fx;fx=pow(1-3/(4*x),1.0/3.0);return(fx);double fnamep(double x) /定义函数 f(x) double fx;fx=pow(1-3/(4*x),-2.0/3.0)*1/(4*pow(x,2.0);return(fx);void main()int n,i;double x,f,fp;x=0.89;/选定初值 n=50;/迭代次数 for(i=1;i=n;i+)/newton迭代 x=x-fname(x)/fnamep(x);coutn=i x=xendl;system(pause);程序清单3.#include#include#include/控制输出有效数字位数函数所在的库using namespace std;void main() static double A10001000;double b1000,x1000,xp1000;int n,i,j,k,m;m=0;double delta,max,eps;delta=0.0;max=1;eps=1e-3;k=0;n=1000;for(i=0;in;i+)/矩阵A赋值 for(j=0;jn;j+)if(i=j)Aij=2.0;elseif(i=j-1)|(i=j+1)Aij=1.0;else Aij=0.0;for(i=0;i=eps)/设置精度 for(i=0;in;i+)/jacobi迭代法 xpi=xi;double ax=0;for(j=0;jn;j+)ax=ax+Aij*xj;ax=ax-Aii*xi;xi=1/Aii*(bi-ax); max=0.0; for(i=0;imax) max=delta; m+;/迭代步数计数 for(i=0;in;i+) coutxi ;coutendl;cout解的精度为setprecision(8)maxendl; max=0.0; for(i=0;imax) max=delta;cout解的误差为setprecision(8)maxendl; cout迭代步数为msyms t f %定义符号变量 f =- 0.000735246*t9 + 0.0359507*t8 - 0.750705*t7 + 8.74608*t6 - 62.32*t5 + 279.926*t4 - 786.973*t3 + 1323.62*t2 - 1192.52*t + 497.288 %定义函数表达式 subs(f,t,18) ans =-4.3447e+06(c)存在龙格现象。对应该列，插值多项式不是一个好的模型，因为当n较大时，该函数在等距节点下的高次插值多项式会发生激烈的震荡。（d）前四组数据插值多项式为： xi=1,2,3,4;yi=67.052,68.008,69.803,72.024;f=lagranz(xi,yi) f =- 0.0688333*t3 + 0.8325*t2 - 1.05967*t + 67.3481998年产量： syms t f %定义符号变量 f=-0.0688333*t3+0.8325*t2-1.05967*t+67.348 %定义函数表达式 subs(f,t,5) ans =74.2580故此时不存在龙格现象。（e）拟合程序： x=1 2 3 4 5 6 7 8 9 10;y=67.052 68.008 69.803 72.024 73.400 72.063 74.669 74.487 74.065 76.777 nafit(x,y,1)ans = 0.9693 66.9034 nafit(x,y,2)ans = -0.0724 1.7652 65.3116 nafit(x,y,3)ans = 0.0182 -0.3729 3.1512 63.7490故一次拟合曲线为：y=0.9693x+66.9034 抛物线拟合曲线为： 三次拟合曲线为： （ii）估计2010年产量：将18带入，得y1=84.3508 y2=73.6276 Y3=105.7934（iii）所以 我觉得 最小二乘拟合直线最能表示表格中的数据。程序清单5.（Simpson）：#include#include#include /控制输出有效数字位数函数所在的库 using namespace std;double fname(double x)/定义函数f(x) double fx;fx=pow(x,3.0)/(pow(x,2.0)+1);return(fx);void main() double a,b,n,h,x,f,eps; cout请输入积分下界a(小数形式)a; cout请输入积分上界b(小数形式)b; cout请输入要求的精度eps(小数形式)eps; couta=a b=b eps=epseps)/控制精度 x=a;f=0.0;/初始化变量 h=(b-a)/(2*n);/步长 for(int i=0;in;i+)/Simpson公式 x=x+i*h; f=f+h/3*(fname(a+2*x)+4*fname(a+2*x+h)+fname(a+2*x+2*h); x=0.0; coutn=n h=setprecision(8)h f=setprecision(12)fendl;n+;/增加节点数 system(pause);程序清单6.（romberg公式）：#include#include#include/控制输出有效数字位数函数所在的库using namespace std;double fname(double x)/定义函数f(x) double fx;fx=pow(x,3.0)/(pow(x,2.0)+1);return(fx);void main()double T100100;double a,b,fn,eps;a=0.0;b=1.0;eps=1e-4;int n,m,i,j;for(i=0;i100;i+) for(j=0;j100;j+) Tij=0;/初始化矩阵 T00=(b-a)/2*(fname(a)+fname(b);n=7;for(j=1;j=n;j+)/计算T0k fn=0.0;for(int p=1;p=(pow(2.0,(double)(j-1);p+)fn=fn+fname(a+(2*p-1)*(b-a)/pow(2.0,(double)j);T0j=T0j-1/2+(b-a)/pow(2.0,(double)j)*fn;for(i=1;i=n;i+)/计算Tml for(j=0;j(n-i);j+)/输出上三角 Tij=(pow(4.0,(double)i)*Ti-1j+1-Ti-1j)/(pow(4.0,(double)i)-1); if(fabs(Ti0-Ti-10)eps)/判断精度 break; for(i=0;in;i+)/输出矩阵 for(j=0;jn;j+) coutsetprecision(12)Tij ; coutendl; system(pause);程序清单七（Eular1)#include#include#include/控制输出有效数字位数函数所在的库using namespace std;double fname(double x,double y) /定义函数f(x，y) double fx;fx=x;return(fx);void main()double x,y,a,b,h;a=0.0;/设置区间 b=1.0;h=0.1; /设置步长 x=a;y=1.0;/给定初值条件 int m,n;m=0;n=(b-a)/h;for(int i=0;i=n;i+)/Euler方法 coutx+i*h setprecision(8)y;y=y+h*fname(x+i*h),y);coutendl;system(pause);程序清单七（Eular2)#include#include#include/控制输出有效数字位数函数所在的库using namespace std;double fname(double x,double y) /定义函数f(x，y) double fx;fx=pow(x,2.0)*y;return(fx);void main() double x,y,a,b,h;a=0.0;/设置区间 b=1.0;h=0.1;/设置步长 x=a;y=1.0;/给定初值条件 int m,n;m=0;n=(b-a)/h;for(int i=0;i=n;i+)/Euler方法 coutx+i*h setprecision(8)y;y=y+h*fname(x+i*h),y);coutendl;system(pause);程序清单七（Eular3)#include#include#include/控制输出有效数