数值分析上机实习报告.docx
- 文档编号:9721053
- 上传时间:2023-02-06
- 格式:DOCX
- 页数:27
- 大小:74.40KB
数值分析上机实习报告.docx
《数值分析上机实习报告.docx》由会员分享,可在线阅读,更多相关《数值分析上机实习报告.docx(27页珍藏版)》请在冰豆网上搜索。
数值分析上机实习报告
(数值分析上机实验报告)
院系:
矿业学院
专业:
矿业工程
班级:
2015
姓名:
王
学号:
2015022
指导教师:
代
第一题
1.用Newton法求解方程X28x4140,在(0.1,1.9)的近似根(初始近似值取为区间端点,迭代6次或误差小于0.00001)。
1.1理论依据及方法应用条件
Newton迭代法:
由一般迭代函数s(x)x
(1)j^)[f(x)][f(x)]j,取s=2时,有2(x)x丄凶j1j!
f(x)
可得二阶迭代序列xk1xk丄型,此种迭代法称为Newton迭代法。
f'(xQ
定理:
设函数在有限区间[a,b]上二阶导数存在,且满足条件
(l)f(a)f(b)0;
(n)f''(x)在区间[a,b]上不变号;
(m)f'(x)工0;
(W)|f(c)|/b-a<|f'(c)|其中c是a,b中使min[|f'(a),f'(b)]达到的一个;
则对任意时近似值x°€[a,b],由Newton迭代过程有:
xk1(xk)xkf(xk)k=0,1,2…
k1kkf'(xk)
所产生的迭代序列{x。
}平方收敛于方程f(x)=0区间[a,b]上的唯一解a。
推论:
设函数f(x)满足定理中条件I,n,m,若选初值x0,使f(x0)•f''(xo)>0,则Newton迭代过程
f(Xk)(k=0,1,2…)产生的迭代序列{Xk}单调收敛于f(x)=0的唯一解a
f'(Xk)
1.2计算程序
#inelude
#inelude
#inelude
#inelude
//牛顿迭代函数
〃牛顿迭代子函数
〃初始数据
usingnamespaeestd;
double*newton(doublea,doubleb,doubleeps);doublenewtonz(doublex);
voidmain()
{
doublea=0.1,b=1.9,eps=0.00001,*result;
eout<<"\n牛顿法解方程:
xA7-28xA4+14=0,在(0.1,1.9)中求近似根,初始值为区间端点,\n误差为0.00001。
\n"< eout<<"学号: 2014021966姓名: 徐林\n"< result=newton(a,b,eps); if(a<=result[0]&&result[0]<=b) cout<<"近似根为: "< if(a<=result[1]&&result[1]<=b) cout<<"近似根为: "< // cout<<"\n"<<"结束,按任意键关闭"< }//主函数结束 //* ****************************************************************** doublenewtonz(doublex){ doublex1=0.0,t; //牛顿迭代子函数 t=(7*pow(x,6)-4*28*pow(x,3)); if(t==0) exit(0); x1=(x-((pow(x,7)-28*pow(x,4)+14)/t)); returnx1; } double*newton(doublea,doubleb,doubleeps){ doublex0=0.0,x1=1.0,x2=0.0,re[2];intk=0; //牛顿迭代函数 x0=a; while(x0>eps) { k++; x2=x1; x1=newtonz(x1);x0=fabs(x1-x2); }re[0]=x1; //代入a迭代计算 //调用牛顿迭代子函数 x0=b,k=0; while(x0>eps) { k++; x2=x1; x1=newtonz(x1); x0=fabs(x1-x2); }re[1]=x1; //代入b迭代计算 //调用牛顿迭代子函数 returnre; 1.3计算结果打印 ■1'C: \U5er5\AdminiStratoDesl'lop、.D亡bug\胡亚梅eye" li.845497 Pi'essanj/ke^tocontinue 1.4MATLAB上机程序 functiony=Newton(f,df,xO,eps,M)d=0; fork=1: M iffeval(df,xO)==O d=2;break else x1=xO-feval(f,xO)/feval(df,xO); end e=abs(x1-x0); x0=x1; ife<=eps&&abs(feval(f,x1))v=epsd=1;break end end ifd==1 y=x1; elseifd==0 y='迭代M次失败: else y='奇异’ end functiony=df(x) y=7*xA6-28*4*xA3; End functiony=f(x) y=xA7-28*xA4+14; End >>x0=1.9; >>eps=0.00001; >>M=100; >>x=Newton('f,'df,x0,eps,M); >>vpa(x,7) 1.5问题讨论 1.需注意的是,要使用Newton迭代法须f(x)x728x414满足定理中的条件I,n,川,以及 f(xc)•f''(x°)>0。 要用误差范围来控制循环的次数,保证循环的次数和质量,编写程序过程中要 注意标点符号的使用,正确运用适当的标点符号,Newton迭代法是局部收敛的,在使用时应先确定 初始值,否则所得的解可能不在所要求的范围内 (3)因为newton法求方程是平方收敛的,所以较为精确,但是要求出函数的导数,且必须有二阶导数。 第二题 2.已知函数值如下表: X 1 2 3 4 5 f(x) 0 0.69314718: 1.0986123 1.3862944 1.6094378 X 6 7 8 9 10 f(x) 1.7917595 1.9459101 2.079445 2.1972246 2.3025851 f'(x) f' (1)=1 f'(10)=0.1 试用三次样条插值求f(4.563)及f'(4.563)的近似值 2.1理论依据及方法应用条件三次样条插值函数可定义为: 对于[a,b]上的一个划分n a 如果定义在[a,b]上的函数S(x),满足 (1).在[xi,xi+i]上为3次多项式; (2).S(x),S'(x), S"(x)在[a,b]上连续,则称S(x)在[a,b]上划分的3次样条函数,如果对于f(x)[a,b], s(x)f(x)还满足s(xjf(x),i0~n,则称s(x)为f(x)的三次样条插值函数。 其基本思想是对均匀分划的插值函数的构造,三次样条函数空间中1,x,,x2,x3,(x-xj)+3为基函数,而取B样条函数Q3((x-Xj)/h)为基函数.由于三次样条函数空间是N+3隹的,故我们把分点扩大到X1,Xn+1则任意三次样条函数可用Q3((X-Xj)/h)线性组合来表示S(x)=N1CjQ3((x-xj)/h)这样对不同插值 j1 问题,若能确定Cj由解的唯一性就能求得S(x)。 由s(xi)=yi,I=1,2,…Ns'(x0)=y0,s(xN)=yn可得 N1 S(Xi)=CjQ3((Xi-Xj)/h)=yi (x°)=1/h CjQ3 (X0-Xj) /h)=y (xN)=1/h CjQ3 (XN-Xj) /h)=y d0 d1 dN d0 A—y'o) 00 d0 yy. (j1jh•4h.( j1j h. j (yN hN1 yN yj hj1 yj1 )6f(Xj1,Xj,Xj1) /*宏定义*/ /*追赶法求解三弯矩方程*/ hj1 hj 2.2计算程序 #include #include #defineN10 main() { floats,ds,t; floatdy0=1,dy9=0.1; intj; intx[N]={1,2,3,4,5,6,7,8,9,10}; floaty[N]={0,0.69314718,1.0986123,1.3862944,1.6094378, 1.7917595,1.9459101,2.079445,2.1972246,2.3025851}; intb[N]={2,2,2,2,2,2,2,2,2,2},h[N-1]; floatd[N],u[N-1],v[N-1],a[N-1],c[N-1],B[N],l[N-1],p[N],X[N]; for(j=1;jv=9;j++) h[j-1]=x[j]-x[j-1]; d[0]=6/h[0]*(y[1]/h[0]-y[0]/h[0]-dy0); d[9]=6/h[8]*(dy9-y[9]/h[8]+y[8]/h[8]); for(j=1;j<=8;j++)d[j]=6/(h[j-1]+h[j])*(y[j+1]/h[j]-y[j]/h[j]-y[j]/h[j-1]+y[j-1]/h[j-1]); for(u[8]=1,j=0;j<=7;j++) u[j]=h[j-1]/(h[j-1]+h[j]); for(v[0]=1,j=1;jv=8;j++) v[j]=h[j]/(h[j-1]+h[j]); for(j=0;j<=8;j++) a[j]=u[j]; for(j=0;j<=8;j++) c[j]=v[j]; for(B[0]=b[0],j=1;jv=9;j++) B[j]=b[j]-a[j]/B[j-1]*c[j-1]; for(j=1;j<=9;j++) l[j]=a[j]/B[j-1]; for(j=1;j<=9;j++) p[j]=d[j]-l[j]*p[j-1]; X[9]=p[9]/B[9]; for(j=8;j>=0;j--) X[j]=p[j]/B[j]-c[j]*X[j+1]/B[j]; t=4.563; /*解f(x)的值 /*解f'(x) /*打印结果*/ s=X[3]*pow((x[4]-t),3)/6/h[3]+X[4]*pow((t-x[3]),3)/6/h[3]+(y[3]-X[3]*h[3]*h[3]/6)*(x[4]/h[3]-t/h[3])+(y[4]-X[4]*h[3]*h[3]/6)*(t/h[3]-x[3]/h[3]); */ ds=-X[3]*pow((x[4]-t),2)/2/h[3]+X[4]*pow((t-x[3]),2)/2/h[3]-(y[3]-X[3]*h[3]*h[3]/6)/h[3]+(y[4]-X[4]*h[3]*h[3]/6)/h[3];的值*/ printf("s=%f\nds=%f\n",s,ds); }2.3计算结果打印 2.4MATLAB上机程序 functionQ=san(ssss,p) Q=zeros(2,1); x=[1;2;3;4;5;6;7;8;9;10]; y=[0;0.69314718;1.0986123;1.3862944;1.6094378;1.7917595;1.9459101;2.079445;2.1972246;2.3025851];h=zeros(10,1); d=zeros(10,1); u=zeros(10,1); v=zeros(10,1); r=zeros(10,1); l=zeros(10,1); z=zeros(10,1); m=zeros(10,1); fort=1: 1: 9; h(t)=x(t+1)-x(t); end d (1)=6/h (1)*((y (2)-y (1))/h (1)-1); d(10)=6/h(9)*(0.1-(y(10)-y(9))/h(9)); fort=1: 1: 8 u(t+1)=h(t)/(h(t)+h(t+1)); v(t+1)=1-u(t+1); d(t+1)=6/(h(t)+h(t+1))*((y(t+2)-y(t+1))/(x(t+2)-x(t+1))-(y(t+1)-y(t))/(x(t+1)-x(t))); end u(10)=1;v (1)=1;r (1)=d (1); fort=2: 1: 10 l(t)=u(t)/r(t-1); r(t)=d(t)-l(t)*v(t-1); end z (1)=d (1); fort=2: 1: 10 z(t)=d(t)-l(t)*z(t-1); end m(10)=z(10)/r(10); fort=9: -1: 1 m(t)=(z(t)-v(t)*m(t+1))/r(t); end fort=1: 1: 10 ifp>=t&&p<(t+1) Q(: 1)=feval(ssss,p,t,x,m,h,y);break end end functionQ=ssss(p,t,x,m,h,y) Q=zeros(2,1); Q(1,1)=((power((x(t+1)-p),3)*m(t)+power((p-x(t)),3)*m(t+1))/6+(y(t)-m(t)*h(t)*h(t)/6)*(x(t+1)-p)+(y(t+1)-m(t+1)*h(t)*h(t)/6)*(p-x(t)))/h(t); Q(2,1)=(-(power((x(t+1)-p),2)*m(t)+power((p-x(t)),2)*m(t+1))/2+(y(t)-m(t)*h(t)*h(t)/6)+(y(t+1)-m(t+1)*h(t)*h(t)/6))/h(t);end 2.5问题讨论 1.若要用追赶法求解三对角方程组,三对角阵需要满足: A(i=1,2,…,n)均非奇异,保证A有唯一的 Doolittle分解;aG丰0; 2.样条插值效果比Lagrange插值好,三次样条插值的解存在且唯一,近似误差较小.并且没有Runge现象。 第二题 3.用Romberg法求: 3 xx1.4(5x7)sinx2dx(允许误差& 3.1理论依据及方法应用条件弹b? f(a)f(b)数值积分的Romberg算法计算步骤如下: ⑴1(I1) T1T1 2 4mTmk)Tm =0.00001)。 (0)(0) 111 时, (k1) Tm1 就停机 ba (k1) 211 fa i1 (2i m1,2,,lk1,2,,lm1 3.2计算程序 /*定义函数f(x)*/ #include { floaty; y=pow(3,x)*pow(x,1.4)*(5*x+7)*sin(x*x); return(y); } main() { floatT[N+1][N+1],h[N+1],a=1,b=3,m[N+1]; inti,l; T[1][0]=(b-a)*(f(a)+f(b))/2; l=1; while(l<=N) { m[l]=0; for(i=1;i<=(pow(2,l-1));i++) m[l]+=f(a+(2*i-1)*(b-a)/pow(2,l)); T[1][l]=(T[1][l-1]+(b-a)*m[l]/pow(2,l-1))/2; l++; } i=1; while(i<=N) { for(l=1;l<=N-i+1;l++) T[i+1][l-1]=(pow(4,i)*T[i][l]-T[i][l-1])/(pow(4,i)-1); h[i]=T[i][0]-T[i+1][0]; if(fabs(h[i])<=1e-5)break; i++; } printf("Theansweris: %f",T[i+1][0]); } 3.3计算结果打印 3.4MATLAB上机程序 function[T,n]=mromb(f,a,b,eps) ifnarginv4,eps=1e-6;end h=b-a; R(1,1)=(h/2)*(feval(f,a)+feval(f,b)); n=1;J=0;err=1; while(err>eps) J=J+1;h=h/2;S=0; fori=1: n x=a+h*(2*i-1); S=S+feval(f,x); end R(J+1,1)=R(J,1)/2+h*S; fork=1: J R(J+1,k+1)=(4Ak*R(J+1,k)-R(J,k))/(4Ak-1); end err=abs(R(J+1,J+1)-R(J+1,J)); n=2*n; end R; T=R(J+1,J+1); formatlong f=@(x)(3.Ax)*(x.A1.4)*(5*x+7)*sin(x*x); [T,n]=mromb(f,1,3,1.e-5) 3.5问题讨论 1、Romberge算法的优点是: a、把积分化为代数运算,而实际上只需求「⑴,以后用递推可得 b、算法简单且收敛速度快,一般4或5次即能达到要求。 c、节省存储量,算出的可存入。 2、Romberge算法的缺点是: a、对函数的光滑性要求较高。 b、计算新分点的值时,这些数值的个数成倍增加。 第四题 4.用定步长四阶Runge-Kutta法求解 厂dy,/dt1 dyz/dty3 dy3/dt10001000y2100y3 yi(0)0 y2(0)0 Iy3(0)0 h0.0005,打印yi(0.025),y(0.045),%(0.085),%(0.1),(i1,2,3) 4.1理论依据及方法应用条件 Runge-Kutta法的基本思想: 旳1不是按Taylor公式展开,而是先写成tn处附近的值的线性组合(有待定系数),再按Taylor公式展开,然后确定待定常数。 四阶古典Runge-Kutta公式: Yn1 Yn(K12K 6 22K3K4) K1 hF(Xn,Yn) K2 hF(Xn £h,Yn 2 2K1) 2 K3 hF(Xn ^h,Yn 2 丄心) 2 K4 hF(Xn £h,Yn 2 1 -K3) 2 4.2计算程序 #include { inti; doubleh=0.0005; doublek1,k2,k3,k4; doubley1=0.0,y2=0.0,y3=0.0; for(i=1;i<=200;i++) { k1=k2=k3=k4=h*1.0; y1+=(k1+2*k2+2*k3+k4)/6; k仁k2=k3=k4=h*y3;y2+=(k1+2*k2+2*k3+k4)/6; k仁h*(1000-1000*y2-100*y3);k2=h*(1000-1000*y2-100*(y3+0.5*k1));k3=h*(1000-1000*y2-100*(y3+0.5*k2));k4=h*(1000-1000*y2-100*(y3+k3));y3+=(k1+2*k2+2*k3+k4)/6; if(i==50) { printf("\ny1(0.025)=%fcontinue; } if(i==90) { printf("\ny1(0.045)=%fcontinue; } if(i==170) { printf("\ny1(0.085)=%fcontinue; } if(i==200) printf("\ny1(0.100)=%f } y2(0.025)=%f y2(0.045)=%f y2(0.085)=%f y2(0.100)=%f y3(0.025)=%f",y1,y2,y3); y3(0.045)=%f",y1,y2,y3); y3(0.085)=%f",y1,y2,y3); y3(0.100)=%f\n\n",y1,y2,y3); 4.3计算结果打印 4.4MATLAB上机程序 functionY=R_K(df1,a,b,h) m=(b-a)/h; Y=zeros(3,1); S=zeros(3,1); K=zeros(3,4); x=a;y1=a;y2=a;y3=a; forn=1: m K(: 1)=feval(df1,x,y1,y2,y3); x=x+0.5*h;S(: 1)=Y+0.5*h.*K(: 1); y1=S(1,1);y2=S(2,1);y3=S(3,1); K(: 2)=feval(df1,x,y1,y2,y3); S(: 1)=Y+0.5*h.*K(: 2); y1=S(1,1);y2=S(2,1);y3=S(3,1); K(: 3)=feval(df1,x,y1,y2,y3); x=x+0.5*h;S(: 1)=Y+h.*K(: 3); y1=S(1,1);y2=S(2,1);y3=S(3,1); K(: 4)=feval(df1,x,y1,y2,y3); Y=Y+h.*(K(: 1)+2.*K(: 2)+2.*K(: 3)+K(: 4))/6; end functionZ=df1(x,y1,y2,y3) Z=zeros(3,1); Z (1)=0*x+0*y1+0*y2+0*y3+1; Z (2)=0*x+0*y1+0*y2+1*y3; Z(3)=0*x+0*y1-1000*y2-100*y3+1000; end 4.5问题讨论 1.定步长四阶runge-kutta法稳定,精度高,可根据有y'f(t,y)变化的情况与需要的精度自动修 改步长,误差小且程序简单,存储量少 2.但是Runge-Kutta法需要每步都计算函数值f(t,y)四次,在函数较复杂时,工作量就会变得较大 可靠性有待核查。 第五题 5.已知A与b
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 数值 分析 上机 实习 报告