数值分析上机题.doc

精选优质文档-倾情为你奉上数值分析上机题习题117（上机题）舍入误差与有效数设，其精确值为。（1）编制按从大到小的顺序，计算的通用程序。（2）编制按从小到大的顺序，计算的通用程序。（3）按两种顺序分别计算，并指出有效位数。（编制程序时用单精度）（4）通过本上机题你明白了什么？专心-专注-专业按从大到小的顺序计算的通用程序为：#include<iostream.h>float sum(float N)float j,s,sum=0;for(j=2;j<=N;j+) s=1/(j*j-1);sum+=s;return sum; 按从小到大的顺序计算的通用程序为：#include<iostream.h>float sum(float N)float j,s,sum=0;for(j=N;j>=2;j-)s=1/(j*j-1);sum+=s;return sum;从大到小的顺序的值从小到大的顺序的值精确值有效位数从大到小从小到大0.0.740050.650.0.74990.7499440.0.0.36通过本上机题，看出按两种不同的顺序计算的结果是不相同的，按从大到小的顺序计算的值与精确值有较大的误差，而按从小到大的顺序计算的值与精确值吻合。从大到小的顺序计算得到的结果的有效位数少。计算机在进行数值计算时会出现“大数吃小数”的现象，导致计算结果的精度有所降低，我们在计算机中进行同号数的加法时，采用绝对值较小者先加的算法，其结果的相对误差较小。习题220（上机题）Newton迭代法（1）给定初值及容许误差，编制Newton法解方程根的通用程序。（2）给定方程，易知其有三个根，。1由Newton方法的局部收敛性可知存在，当时，Newton迭代序列收敛于根。试确定尽可能大的。2试取若干初始值，观察当，时Newton序列是否收敛以及收敛于哪一个根。（3）通过本上机题，你明白了什么？解：（1）编制的通用程序：#include<iostream.h>#include<math.h>#define eps 0. /给定容许误差float f(float x) /定义函数f(x)float f;f=x*x*x/3-x; /f(x)的表达式;return(f);float df(float x) /定义函数df(x)，计算f(x)的导函数float df;df=x*x-1; /f(x)导函数的表达式;return (df);void main(void)float x0,x1,a;int k=0;cout<<"请输入初值x0:"cin>>x0;doa=-f(x0)/df(x0);x1=x0+a;k+;x0=x1;while(fabs(a)>eps);cout<<k<<'t'<<x0; /输出迭代的次数和根值（2）计算迭代序列收敛于根的尽可能大的的函数为：#include<iostream.h>#include<math.h>void delay(int n) /定义延时函数for(n=10000;n>0;n-);#define eps 0.float f(float x) /定义函数f(x)float f;f=x*x*x/3-x; /f(x)的表达式;return(f);float df(float x) /定义函数df(x)，计算f(x)的导函数float df;df=x*x-1; /f(x)导函数的表达式;return (df);int judgement(float z)int count=5;float x0,x1,type,type1;x0=z;while(count->0) x1=x0-f(x0)/df(x0);type=fabs(x1); type1=fabs(x1-x0); /调试值用cout<<"count="<<count<<'t'<<"type="<<type<<'t'<<"type1="<<type1<<'n'if(fabs(x1-x0)<eps)return 1;x0=x1;delay(30000); /调试值用 return 0; void main(void)float delta=0;int flag=1;while(flag=1)cout<<"方程的根为:"<<'n'delta+=eps;flag=judgement(delta);cout<<"输出方程根收敛的区间值:n"cout<<delta-eps; /输出收敛的区间值当步长为0.001时，程序计算出的的为=0.774，即在区间（-0.774，0.774）内迭代序列收敛于0。对于不同得初始值收敛于不同的根, 在（-，-1）内收敛于，在（-0.774，0.774）内收敛于，在（1，+）内收敛于，但在内（0.774，1）和（1，0.774）均可能收敛于和。，分别为方程的精确解。分析：对于不同的初值，迭代序列会收敛于不同的根，所以在某个区间内求根对于初值的选取有很大的关系。产生上述结果的原因是区间不满足大范围收敛的条件。习题339 （上机题）列主元Gauss消去法 对于某电路的分析，归结为求解线性方程组RI=V。 其中,31 -13 0 0 0 -10 0 0 0 -13 35 -9 0 -11 0 0 0 00 -9 31 -10 0 0 0 0 00 0 -10 79 -30 0 0 0 -90 0 0 -30 57 -7 0 -5 00 0 0 0 -7 47 -30 0 00 0 0 0 0 -30 41 0 00 0 0 0 -5 0 0 27 -20 0 0 -9 0 0 0 -2 29 R= VT=(-15,27, -23,0,-20,12,-7,7,10)T（1） 编制解n阶线性方程组Ax=b的列主元三角分解法的通用程序；（2）用所编制的程序解线性方程组RI=V，并打印出解向量，保留五位有效数；（3）本编程之中，你提高了哪些编程能力？程序为：#include<iostream.h>#include<math.h>void main(void)int i,j,n,k,q;float a1011,s10,s110;cout<<"请输入n的值:"cin>>n;cout<<"输入数组a:"<<endl;for(i=1;i<=n;i+)for(j=1;j<=(n+1);j+)cin>>aij; /给矩阵a赋值for(i=1;i<=n;i+)for(j=1;j<=(n+1);j+)cout<<aij<<'t'cout<<'n' /输出数组acout<<"'''''''''''''''''''''''''"<<'n'/进行第一行和第一列元素的求取'''''''''''''''''''''''''/int t=1;for(i=1;i<=n;i+)si=ai1;float max=fabs(s1);for(i=2;i<=n;i+)if(fabs(si)>max)max=fabs(si);t=i;for(j=1;j<=(n+1);j+)float b=a1j;a1j=atj;atj=b; /进行第一列主元互换for(i=2;i<=n;i+)ai1=ai1/max; /第一列除以a11for(i=1;i<=n;i+)for(j=1;j<=(n+1);j+)cout<<aij<<'t'cout<<'n'/输出进行第一步变换的数组acout<<"'''''''''''''''''''''''''"<<'n'/进行第k步分解'''''''''''''''''''''''''''''''''''''''''/for(k=2;k<=n;k+)for(i=k;i<=n;i+)float sum=0;for(q=1;q<k;q+)sum+=aiq*aqk;s1i=aik-sum;int l=k;float m=fabs(s1k);for(i=k;i<=n;i+)/比较第k步分解的第k列值的大小if(fabs(s1i)>m)m=fabs(s1i);l=i; /返回行值for(j=1;j<=n+1;j+) /交换两行元素float s2=akj;akj=alj;alj=s2;for(j=k;j<=n+1;j+)/算出第k行行元素的值float sum1=0;for(q=1;q<k;q+)sum1+=akq*aqj;akj=akj-sum1;for(i=k+1;i<=n;i+)/算出第k列列元素的值float sum2=0;for(q=1;q<k;q+)sum2+=aiq*aqk;aik=(aik-sum2)/(akk); /第k步分解结束for(i=1;i<=n;i+)for(j=1;j<=(n+1);j+)cout<<aij<<'t'cout<<'n' /输出改变后的数组/输出解'''''''''''''''''''''''''''''''''''''''''''''''''''''''/float x10;for(i=n-1;i>=1;i-)xn=ann+1/ann;float sum3=0;for(j=i+1;j<=n;j+)sum3+=aij*xj;xi=(ain+1-sum3)/aii; /回代过程for(i=1;i<=n;i+)cout<<'x'<<i<<'='<<xi<<endl; /输出解向量结果：方程的解为：x1= -0.28923，x2= 0.34544，x3= -0.71281，x4= -0.22061，x5= -0.43040，x6= 0.15431，x7= -0.，x8= 0.20105，x9= 0.29023。分析：我感觉是提高了查错误点的能力和编循环语句的能力，即利用很规整的迭代公式进行编程。另外列主元三角分解法的阶梯步骤有了更深的了解!习题 437.（上机题）3次样条插值函数（1）编制求第一型3次样条插值函数的通用程序; （2） 已知汽车曲线型值点的数据如下：0123456789102.513.304.044.705.225.545.785.405.575.705.80端点条件为=0.8，=0.2。用所编制程序求车门的3次样条插值函数S(x),并打印出S(i+0.5)(i=0,1,9)。解：通用程序：#include<iostream.h>void main(void)float x11;/存放数组xjfloat y11;/存放数组yjfloat h11;/存放数组hjfloat u11;/存放数组ujfloat v11;/存放数组vjfloat d11;/存放数组djfloat M11;/存放数组Mjfloat b11;/ 存放数组bjfloat t11,l11,yy11,s4,aa1,aa2,aa3,aa4;float s110;int i,j,n;float xx;/x为区间值/将初值初始化cout<<"请输入n的值:n"cin>>n;cout<<"输入数组x:n"for(i=0;i<=n;i+)cin>>xi;cout<<"输入数组y:n"for(i=0;i<=n;i+)cin>>yi;/输入端点值float df2;cout<<"输入两个端点值:n"for(i=0;i<2;i+)cin>>dfi;/利用书本上的算法求出所需要的值/求出hj的值for(j=0;j<=n-1;j+)hj=xj+1-xj;cout<<'h'<<''<<j<<''<<'='<<hj<<'t'cout<<endl;/求出uj和vj的初值v0=1;un=1;for(j=1;j<=n-1;j+)uj=hj-1/(hj-1+hj);vj=hj/(hj-1+hj);/求出dj的值for(j=1;j<n;j+)dj=6*(yj+1-yj)/hj-(yj-yj-1)/hj-1)/(hj+hj-1);d0=6*(y1-y0)/h0-df0)/h0;dn=6*(df1-(yn-yn-1)/hn-1)/hn-1;for(j=1;j<=n;j+)cout<<'u'<<''<<j<<''<<'='<<uj<<'t'cout<<endl;for(j=0;j<n;j+)cout<<'v'<<''<<j<<''<<'='<<vj<<'t'cout<<endl;for(j=0;j<=n;j+)cout<<'d'<<''<<j<<''<<'='<<dj<<'t'cout<<endl;/利用书本上的追赶法求解方程组for(i=0;i<=n;i+)bi=2;cout<<endl;t0=b0;yy0=d0;/消元过程for(i=1;i<=n;i+)li=ui/ti-1;ti=bi-li*vi-1;yyi=di-li*yyi-1;/回代过程Mn=yyn/tn;for(i=n-1;i>=0;i-)Mi=(yyi-vi*Mi+1)/ti;/将Mj的值输出for(i=0;i<=n;i+)cout<<'M'<<''<<i<<''<<'='<<Mi<<endl;/输出插值多项式的系数for(j=0;j<n;j+)s0=yj;s1=(yj+1-yj)/hj-(Mj/3+Mj+1/6)*hj;s2=Mj/2;s3=(Mj+1-Mj)/(6*hj);cout<<"当x的值在区间"<<'x'<<''<<j<<''<<"到"<<'x'<<''<<(j+1)<<''<<"时,输出插值多项式的系数:n"for(int k=0;k<4;k+)cout<<'s'<<''<<k<<''<<'='<<sk<<'t'cout<<endl;（2）编制的程序求车门的3次样条插值函数S(x)：x属于区间0,1时；S(x)=2.51+0.8(x)-0.(x)(x)-0.(x)(x)(x)x属于区间1,2时；S(x)=3.3+0.(x-1)-0.(x-1)(x-1)-0.(x-1)(x-1)(x-1)x属于区间2,3时；S(x)=4.04+0.(x-2)-0.(x-2)(x-2)-0.(x-2)(x-2)(x-2)x属于区间3,4时；S(x)=4.7+0.(x-3)-0.(x-3)(x-3)-0.(x-3)(x-3)(x-3)x属于区间4,5时；S(x)=5.22+0.(x-4)-0.(x-4)(x-4)+0.(x-4)(x-4)(x-4)x属于区间5,6时；S(x)=5.54+0.(x-5)+0.(x-5)(x-5)-0.(x-5)(x-5)(x-5)x属于区间6,7时；S(x)=5.78-0.(x-6)-0.(x-6)(x-6)+0.(x-6)(x-6)(x-6)x属于区间7,8时；S(x)=5.4-0.(x-7)+0.(x-7)(x-7)-0.(x-7)(x-7)(x-7)x属于区间8,9时；S(x)=5.57+0.(x-8)-0.(x-8)(x-8)+0.(x-8)(x-8)(x-8)x属于区间9,10时；S(x)=5.7+0.(x-9)-0.(x-9)(x-9)+0.(x-9)(x-9)(x-9)S(0.5)=2.90856 S(1.5)=3.67843 S (2.5)=4.38147S(3.5)=4.98819 S(4.5)=5.38328 S(5.5)=5.7237S(6.5)=5.59441 S(7.5)=5.42989 S(8.5)=5.65976S(9.5)=5.7323 习题623（上机题）常微分方程初值问题数值解（1）编制RK4方法的通用程序；（2）编制AB4方法的通用程序（由RK4提供初值）；（3）编制AB4-AM4预测校正方法的通用程序（由RK4提供初值）；（4）编制带改进的AB4-AM4预测校正方法的通用程序（由RK4提供初值）;（5）对于初值问题取步长,应用（1）（4）中的四种方法进行计算，并将计算结果和精确解作比较；（6）通过本上机题，你能得到哪些结论？解：程序为：#include<iostream.h>#include<fstream.h>#include<stdlib.h>#include<math.h>ofstream outfile("data.txt");/此处定义函数f(x,y)的表达式/用户可以自己设定所需要求得函数表达式double f1(double x,double y)double f1;f1=(-1)*x*x*y*y;return f1;/此处定义求函数精确解的函数表达式double f2(double x)double f2;f2=3/(1+x*x*x);return f2;/此处为精确求函数解的通用程序void accurate(double a,double b,double h)double x100,accurate100;x0=a;int i=0;outfile<<"输出函数准确值的程序结果:n"doxi=x0+i*h;accuratei=f2(xi);outfile<<"accurate"<<i<<"="<<accuratei<<'n'i+;while(i<(b-a)/h+1);/此处为经典Runge-Kutta公式的通用程序void RK4(double a,double b,double h,double c)int i=0;double k1,k2,k3,k4;double x100,y100;y0=c;x0=a;outfile<<"输出经典Runge-Kutta公式的程序结果:n"doxi=x0+i*h;k1=f1(xi,yi);k2=f1(xi+h/2),(yi+h*k1/2);k3=f1(xi+h/2),(yi+h*k2/2);k4=f1(xi+h),(yi+h*k3);yi+1=yi+h*(k1+2*k2+2*k3+k4)/6;outfile<<"y"<<""<<i<<"="<<yi<<'n'i+;while(i<(b-a)/h+1);/此处为4阶Adams显式方法的通用程序void AB4(double a,double b,double h,double c)double x100,y100,y1100;double k1,k2,k3,k4;y0=c;x0=a;outfile<<"输出4阶Adams显式方法的程序结果:n"for(int i=0;i<=2;i+)xi=x0+i*h;k1=f1(xi,yi);k2=f1(xi+h/2),(yi+h*k1/2);k3=f1(xi+h/2),(yi+h*k2/2);k4=f1(xi+h),(yi+h*k3);yi+1=yi+h*(k1+2*k2+2*k3+k4)/6;int j=3;y10=y0;y11=y1;y12=y2;y13=y3;doxj=x0+j*h;y1j+1=y1j+(55*f1(xj,y1j)-59*f1(xj-1,y1j-1)+37*f1(xj-2,y1j-2)-9*f1(xj-3,y1j-3)*h/24;outfile<<"y1"<<""<<j<<"="<<y1j<<'n'j+;while(j<(b-a)/h+1);/主函数void main(void)double a,b,h,c;cout<<"输入上下区间、步长和初始值:n"cin>>a>>b>>h>>c;accurate(a,b,h);RK4(a,b,h,c);AB4(a,b,h,c);结果为：由经典Runge-Kutta公式得出的结果列在下面的表格中，以及精确值y(xi)和精确值和数值解的误差：ixiyiy(xi)|y(xi)-yi|0033010.12.9972.9971.87138e-00720.22.976192.976193.91665e-00730.32.921132.921137.58342e-00740.42.819552.819551.61101e-00650.52.666662.666673.17735e-00660.62.46712.467115.00551e-00670.7 2.23382.23385.77233e-00680.81.984121.984134.12954e-00690.91.735111.735111.15554e-007101.01.500011.55.80668e-006111.11.287011.2871.13075e-005121.21.099721.099711.54242e-005131.30.0.938381.77272e-005141.40.80130.1.83754e-005151.50.0.1.78e-005由AB4方法得出的结果为：Y10=3 y11=2.997 y12=2.97619 y13=2.92113 y14=2.81839y15=2.66467 y16=2.4652 y17=2.23308 y18=1.98495 y19=1.73704y110=1.50219 y111=1.28876 y112=1.10072 y113=0.93871 y114=0.y115=0.通过本上机题我明白了各种求微分方程的数值方法，经典Runge-Kutta公式，AB4方法以及AB4-AM4预测校正方法求解公式的精度是不同的。其中经典Runge-Kutta公式的精度，四阶Adams显式方法（AB4）具有4阶精度。 习题710.抛物线方程Crank-Nicolson格式（1）编制用Crank-Nicolson格式求抛物线方程2u/2x = f(x,t) (0<x<1, 0u(x,0) = (0u(0,t) = ，u（1,t）= (0<t数值解的通用程序。（2）就a=1, f(x,t)=0,=exp(x),=exp(t),=exp(1+t),M=40,N=40,输入点（0.2，1.0），（0.4，1.0），（0.6，1.0），（.8，1.0）4点处u（x,t）的近似值。（3） 已知所给方程的精确解为u(x,t)=exp(x+t),将步长反复二分，从（0.2，1.0），（0.4，1.0），（0.6，1.0），（0.8，1.0）4点处精确解与数值解的误差观察当步 长缩小一半时，误差以什么规律缩小。算法： 本题选择空间步长h=0.025,时间步长=0.025 1. 置初值； 2.通过点(i-1,k),(i-1,k+1),(i,k),(i+1,k),(i+1,k+1)求解点(I,k+1)的值； 3.输出结果 原程序：#include <iostream>#include <math.h>float h=0.025,k=0.025;int m=40;int n=40; float y4040,r=a*k/(h*h);void Input() int i,j; cout<<"Loading Input Data."<<endl; for(i=0;i<m;i+) for(j=0;j<n;j+) if (i=j) aij=1+r; for(j=0;j<n;j+) if (j=i+1)|(i=j+1) aij=-r/2;int main() Input(); /read data int r,i; for(k=0;k<(m-1);k+) int select(); /select main element r=select(); void exchange(int g); exchange(r); /exchange void analyze(); analyze(); /analyze void ret(); ret(); / replace back cout<<"The solution vector is below:"<<endl; for(i=0;i<m;i+) cout<<"x"<<i<<"="<<xi<<endl; return 0;int select() int f,t; float max; f=k; float si(int u,int v); max=float(fabs(si(k,k); for(t=(k+1);t<(m-1);t+) if(max<fabs(si(t,k) max=float(fabs(si(t,k); f=t; return f;float si(int u,int v) float sum=0; int q; for(q=0;q<k;q+) sum+=auq*aqv; sum=auv-sum; return sum;void exchange(int g) int t; float temp; for(t=0;t<n;t+) temp=akt; akt=agt; agt=temp; void analyze() int t; float si(int u,int v); for(t=k;t<n;t+) akt=si(k,t); for(t=(k+1);t<m;t+) atk=(float)(si(t,k)/akk);void ret() int t,z;float sum; xm-1=(float)am-1m/am-1m-1; for(t=(m-2);t>-1;t-) sum=0; for(z=(t+1);z<m;z+) sum+=atz*xz; xt=(float)(atm-sum)/att; 运行结果: =Numerical Solution of Partial Differential Equations= NumericalResult AccurateValue Errors（0.2,1.0） 3.31947 3.32012 0.00045（0.4,1.0） 4.05432 4.05520 0.00088（0.6,1.0） 4.95238 4.95303 0.00065（0.8,1.0） 4.04897 6.04965 0.00067分析与结论: Crank-Nicolson格式，步长越小，精度高。但计算过程比较麻烦。