数值分析上机题参考答案.docx
研究生数值分析上机试题及解答

东华大学研究生数值分析试题(上机部分)A 卷2008年12月 时间:60分钟班级 学号 机号 姓名 得分 注意:要求写出M 函数(如果需要)、MATLAB 命令和计算结果。
1. 求下列方程组在0<α, β<1中的解⎩⎨⎧-=+=βαββααsin 2.0cos 7.0cos 2.0sin 7.0 命令fun=inline('[x(1)-0.7*sin(x(1))-0.2*cos(x(2)),x(2)-0.7*cos(x(1))+0.2*sin(x(2))]','x'); [x,f,h]=fsolve(fun,[0.5 0.5]) 结果α=0.5265,β=0.50792命令>> fun=inline('c(1)+c(2)*x.^2','c','x'); >> x=[1.1 1.3 1.4 1.6 1.8]; >> y=[26 22 23 24 25];>> c=lsqcurvefit(fun,[0 0],x,y) 结果 c =23.7256 0.12873.求解下列微分方程组2(0)2013(0)1x x yx t y x yy '=-=⎧<<⎨'=+=⎩(结果只要求写出t =1时的解) 命令>> fun=inline('[y(1)-2*y(2);3*y(1)+y(2)]','t','y'); >> [t,y]=ode45(fun,[0 1], [2 1]) 结果x(1)=-5.6020, y(1)=2.15634.用定步长Gauss 积分法(课本123页)计算积分31e ln(1)x x dx -+⎰的近似值(等分数取4,每段取2个Gauss 点)。
命令fun=inline('exp(-x).*log(1+x)','x'); nagsint(fun,1,3,4,2) 结果 0.30865.矩阵改进平方根分解(课本25页)的计算公式为: d 1=a 11, 对i =2, 3, ⋯, n ,iki k ik ii i j ij ij j k jk ik ij ij l s a d i j d s l l s a s ∑∑-=-=-=-==-=1111,1,,2,1 ,/ ,试编写矩阵改进平方根分解的程序,并求矩阵1111551514A -⎛⎫ ⎪=-- ⎪ ⎪-⎝⎭的改进平方根分解。
数值分析上机答案computer中科院

源程序如下:function[x_jac,x_gauss,x_sor,A,b]=diedai(n,w,e,N,Ain,bin) if nargin<3e=1e-5;N=1000;endif nargin<4N=1000;endif nargin<2error('²ÎÊý²»¹»');endA=hilb(n);b=A*ones(n,1);if nargin>4if nargin==5error('ÇëÊäÈëb');elseA=Ain;b=bin;endend%%Ñſ˱ȵü´údiedai_n=0;x=zeros(n,2);while max(abs(b-A*x(:,2)))>ediedai_n=diedai_n+1;if diedai_n>Ndisp('Ñſ˱ȵü´úδ´ïµ½¾«¶È');break;endx(:,1)=x(:,2);for i=1:nif i>1x(i,2)=(b(i)-A(i,1:(i-1))*x(1:(i-1),1)-A(i,(i+1):n)*x((i+ 1):n,1))/A(i,i);elsex(i,2)=(b(i)-A(i,(i+1):n)*x((i+1):n,1))/A(i,i);endendendx_jac.x=x(:,2);x_jac.n=diedai_n-1;x_jac.error=max(abs(b-A*x(:,2)));%%¸ß˹µü´údiedai_n=0;x=zeros(n,2);while max(abs(b-A*x(:,2)))>ex(:,1)=x(:,2);diedai_n=diedai_n+1;if diedai_n>Ndisp('¸ß˹¡ªÈüµÂ¶ûµü´úδ´ïµ½¾«¶È');break;endfor i=1:nif i>1x(i,2)=(b(i)-A(i,1:(i-1))*x(1:(i-1),2)-A(i,(i+1):n)*x((i+ 1):n,1))/A(i,i);elsex(i,2)=(b(i)-A(i,(i+1):n)*x((i+1):n,1))/A(i,i); endendendx_gauss.x=x(:,2);x_gauss.n=diedai_n-1;x_gauss.error=max(abs(b-A*x(:,2)));%%SORµü´údiedai_n=0;x=zeros(n,2);while max(abs(b-A*x(:,2)))>ex(:,1)=x(:,2);diedai_n=diedai_n+1;if diedai_n>Ndisp('SORµü´úδ´ïµ½¾«¶È');break;endfor i=1:nif i>1x(i,2)=x(i,1)+w*(b(i)-A(i,1:(i-1))*x(1:(i-1),2)-A(i,i:n)* x(i:n,1))/A(i,i);elsex(i,2)=x(i,1)+w*(b(i)-A(i,i:n)*x(i:n,1))/A(i,i);endendendx_sor.x=x(:,2);x_sor.n=diedai_n-1;x_sor.error=max(abs(b-A*x(:,2)));end运行结果:由于SOR的松弛因子取1时与高斯-赛德尔迭代相同,故在矩阵n取8,10时省去了SOR迭代松弛因子取1的步骤。
数值分析上机试题及解答2008

东华大学数值分析试题(上机部分)A 卷2008年12月 时间:60分钟班级 学号 机号 姓名 得分 注意:要求写出M 函数(如果需要)、MATLAB 命令和计算结果。
1. 求下列方程组在0<α, β<1中的解⎩⎨⎧-=+=βαββααsin 2.0cos 7.0cos 2.0sin 7.0 命令fun=inline('[x(1)-0.7*sin(x(1))-0.2*cos(x(2)),x(2)-0.7*cos(x(1))+0.2*sin(x(2))]','x'); [x,f,h]=fsolve(fun,[0.5 0.5]) 结果α=0.5265,β=0.50792命令>> fun=inline('c(1)+c(2)*x.^2','c','x'); >> x=[1.1 1.3 1.4 1.6 1.8]; >> y=[26 22 23 24 25];>> c=lsqcurvefit(fun,[0 0],x,y) 结果 c =23.7256 0.12873.求解下列微分方程组2(0)2013(0)1x x yx t y x yy '=-=⎧<<⎨'=+=⎩(结果只要求写出t =1时的解) 命令>> fun=inline('[y(1)-2*y(2);3*y(1)+y(2)]','t','y'); >> [t,y]=ode45(fun,[0 1], [2 1]) 结果x(1)=-5.6020, y(1)=2.15634.用定步长Gauss 积分法(课本123页)计算积分31e ln(1)x x dx -+⎰的近似值(等分数取4,每段取2个Gauss 点)。
命令fun=inline('exp(-x).*log(1+x)','x'); nagsint(fun,1,3,4,2) 结果 0.30865.矩阵改进平方根分解(课本25页)的计算公式为: d 1=a 11, 对i =2, 3, ⋯, n ,iki k ik ii i j ij ij j k jk ik ij ij l s a d i j d s l l s a s ∑∑-=-=-=-==-=1111,1,,2,1 ,/ ,试编写矩阵改进平方根分解的程序,并求矩阵1111551514A -⎛⎫ ⎪=-- ⎪ ⎪-⎝⎭的改进平方根分解。
数值计算方法上机实验考试答案

数值计算方法上机实验考试答案1. 小型火箭初始质量为900千克,其中包括600千克燃料。
火箭竖直向上发射时燃料以15千克/秒的速率燃烧掉,由此产生30000牛顿的恒定推力。
当燃料用尽时引擎关闭。
设火箭上升的整个过程中,空气阻力与速度平方成正比,比例系数为0.4(千克/米)。
重力加速度取9.8米/秒2.A. 建立火箭升空过程的数学模型(微分方程);B. 求引擎关闭瞬间火箭的高度、速度、加速度,及火箭到达最高点的时间和高度。
解:A . 建立模型:)(t x ——t 时刻的火箭高度;T =30000(牛顿)——火箭推力,当时间t >40秒时,T=0; m t m 150-=——火箭飞行过程中的质量,t >40秒时,300=m 千克0m =900(千克)——火箭初始质量; 1m =600(千克)——燃料质量;c =15(千克/秒)——燃料的燃烧速率; k =0.4(千克/米)——空气阻力系数; g =9.8(米/秒2)——重力加速度由能量守恒定律,可得到火箭飞行过程的方程:mg T t x k t x m -+'-=''2)]([)(这是一个初值问题,初始条件为 0)0(,0)0(='=x x设)(),(21t x x t x x '==,则问题化为求下列微分方程组的初值问题:⎪⎩⎪⎨⎧--+--='='g t m T x t m k x x x 151******** 0)0(,0)0(21==x xB . 关闭引擎时4015/600==t (秒),所求的是此时火箭的高度)40(1x ;速度)40(2x ;加速度)40()40(21x x '''或,及火箭到最高点的时间m t 和高度)(1m t x 。
具体的Matlab 程序如下: 首先建立微分方程的的m 文件:function y=huojian(t,x)k=0.4;g=9.8;m0=900;T=30000;m=m0-15*t;if t>40T=0;m=300;endy=[x(2),-(k/m)*x(2)^2+T/m-g]';主程序:%feixing.mk=0.4;g=9.8;m0=900;T=30000;x0=[0,0];ts=0:1:55;[t,x]=ode23('huojian',ts,x0);[t,x(:,1)]%------a=[t,x];x40=a(41,2) %燃料用尽时的高度v40=a(41,3) %燃料用尽时的速度a40=-(k/300)*v40^2+T/300-g %燃料用尽时的加速度%-------xmax=max(x(:,1)) %火箭到达最高点的高度subplot(2,1,1),plot(t,x(:,1)),title('altitude') subplot(2,1,2),plot(t,x(:,2)),title('speed')运行结果为:1.0e+003 *0 00.0010 0.01180.0020 0.04750.0030 0.10670.0040 0.18890.0050 0.29270.0060 0.41680.0070 0.55920.0080 0.71800.0090 0.89140.0100 1.07700.0380 7.80510.0390 8.06310.0400 8.32240.0410 8.54000.0420 8.70040.0430 8.82420.0440 8.92180.0450 8.99940.0460 9.06070.0470 9.10830.0480 9.14370.0490 9.16810.0500 9.18220.0510 9.18640.0520 9.18070.0530 9.16510.0540 9.13900.0550 9.1018x40 =8322.4v40 = 254.1728a40 = 4.0616xmax =9186.4关闭引擎时4015/600==t (秒),此时火箭的高度h=)40(1x =8322.4米,速度v=)40(2x =254.1728米/秒,加速度为a=)40(1x ''= 4.0616米/秒2,火箭到最高点的时间m t =51秒,高度)(1m t x =9186.4米。
数值分析韩旭里答案

数值分析韩旭里答案【篇一:数值分析上机题目】>1631110xxxx 材料科学与工程学院一.第2章插值法l2.7 给定数据表2-15.用newton插值公式计算3次插值多项式n3(x).表2-15x f(x)1 1.251.52.500 1.002 5.50a. matlab代码如下,two.m,%第二章,p45,练习题2第七题 clear(); x=[1,1.5,0,2];y(:,1)=[1.25,2.50,1.00,5.50];%已知点集合x和y syms t w; w(1)=1; %计算基函数序列w和差商表y,以及函数序列的权数diag(y),计算的牛顿三次多项式表述为t的函数 for j=2:length(x) fori=j:length(x)y(i,j)=(y(i,j-1)-y(i-1,j-1))/(x(i)-x(i-j+1)); i=i+1; endw(j)=prod(t-x(1:j-1)); j=j+1; enddisp(三次牛顿插值多项式为); disp(collect(w*diag(y)));plot(x,y(:,1),*); hold on;fplot(collect(w*diag(y)),[-0.5,2.5]);legend({已知点集,三次牛顿插值多项式函数},location,northwest,fontsize,14); xlabel(x,fontsize,16);ylabel(y,fontsize,16); hold off;b. 计算结果如下:二.第3章函数逼近与数据拟合a. matlab代码,three.m,%第三章函数逼近与数据拟合,p68练习题,第2题 clear(); syms x;%所使用的非线性基函数序列,用符号表示 y=abs(x);%被逼近函数f=[1,x^2,x^4];%求解法方程的系数矩阵a*gn=b,其中a和b均为行向量gn=ones(length(f),length(f)); for i=1:length(f) for j=1:length(f) gn(i,j)=int(f(i)*f(j),-1,1);j=j+1; endb(i)=int(f(i)*y,-1,1); i=i+1; enda=b/gn;%最佳平方逼近的系数行向量 disp(逼近函数表达式);disp(vpa(f*a));disp(最佳函数逼近得平方误差); disp(vpa(int(y^2,-1,1)-a*b));fplot(y,[-1,1]); hold on; fplot(a*f,[-1,1]);legend({被逼近函数,逼近函数},location,north,orientation,horizontal,fontsize,16,fontweight,b old);xlabel(x,fontsize,20,fontweight,bold);ylabel(y,fontsize,20,fontweight,bold); hold off;b. 运行结果如下:三.第4章数值积分与数值微分例4.9用romberg求积法计算定积分 01sin?(??)??a. matlab代码,four.m%romberg求积公式,外推原理 clear(); clear(); format long; a=0; b=1;t(1,1)=(b-a)/2*(f(a)+f(b));t(2,1)=1/2*t(1,1)+(b-a)/2*f((a+b)/2); t(1,2)=(4*t(2,1)-t(1,1))/(4-1);col=2;while abs(t(1,col)-t(1,col-1))0.5*10^-6%t(1,col)对应的计算的是多少步的值,col→coln关系col=col+1;%此时求得是第n+1次均分后的结果,使用的是第n次的结果,注意在矩阵 %计算的第n斜列是第n-1次均分的结果 for j=1:colif j==1h=(b-a)/2^(col-2);%使用n+1之前的第n次结果【篇二:数值分析a教学】>一、课程基本信息二、课程目的和任务“数值分析”是理工科院校计算数学、力学、物理、计算机软件等专业的学生必须掌握的一门重要的基础课程。
《数值分析》上机习题

《数值分析》上机习题1.用Newton法求方程X7 - 28X4 + 14=0在(0.1,1.9)中的近似根(初始近似值取为区间端点,迭代6次或误差小于0.00001)。
#include<stdio.h>#include<math.h>int main(){ float x1,x,f1,f2;static int count=0;x1=0.1 ;//定义初始值do{x=x1;f1=x*(x*x*x*x*x*x-28*x*x*x)+14;f2=7*x*x*x*x*x*x-4*28*x*x*x;//对函数f1求导x1=x-f1/f2; count++; }while(fabs(x1-x)>=1e-5);printf("%8.7f\n",x1); printf("%d\n",count);return 0;}2.已知函数值如下表:试用三次样条插值求f(4.563)及f’(4.563)的近似值。
include <stdio.h>#include <math.h>#define N 11main(){double B[N+1][N+1],m,x,u[N+1],y[N+1],c[N+1],d[N+1];double e[N+1]={2,0,4.15888308,6.5916738,8.3177664,9.6566268,10.750557,11.675460 6,12.47667,13.1833476,13.8155106,14.0155106};int i;x=4.563;B[0][0]=-2;B[0][1]=-4;B[N][N-1]=4;B[N][N]=2;for(i=1;i<N;i++){B[i][i-1]=1;B[i][i]=4;B[i][i+1]=1;}u[0]=B[0][0];y[0]=e[0];for(i=1;i<N;i++){m=B[i][i-1]/u[i-1];u[i]=B[i][i]-m*B[i-1][i];y[i]=e[i]-m*y[i-1];}c[N]=y[N]/u[N];for(i=N-1;i>=0;i--)c[i]=(y[i]-B[i][i+1]*c[i+1])/u[i];for(i=0;i<12;i++){m=fabs(x-i);if(m>=2)d[i]=0;else if(m<=1)d[i]=0.5*fabs(pow(m,3))-m*m+2.0/3;elsed[i]=(-1.0/6.0)*fabs(pow(m,3))+m*m-2*fabs(m)+4/3.0;} m=0;for(i=0;i<12;i++) m=m+c[i]*d[i];printf("f(%4.3f)=%f\n",x,m);printf("f'(4.563)=%lf\n",(c[4]-c[2])/2); }3.用Romberg 算法求)00001.0(sin )75(32314.1=+⎰ε允许误差dx x x x x . #include "stdafx.h" #include<stdio.h> #include<math.h> float f(float x) {float f=0.0;f=pow(3.0,x)*pow(x,1.4)*(5*x+7)*sin(x*x); return (f); } main() {int i=1,j,k,n=12;float T[12],a=1.0,b=3.0,s=0.0; T[0]=0.5*(b-a)*(f(a)+f(b));for(j=1;j<n-1;j++){ for(k=1;k<=pow(2,j-1);k++) s+=f(a+(2*k-1)*(b-a)/pow(2,j)); T[j]=0.5*(T[j-1]+(b-a)*s/pow(2,j-1)); s=0.0; }T[11]=(4*T[1]-T[0])/(float)3;for(;fabs(T[11]-T[0])>0.00001;i++) {T[0]=T[11];for(j=1;j<n-1-i;j++) T[j]=(pow(4,i)*T[j+1]-T[j])/(pow(4,i)-1); T[11]=(pow(4,i+1)*T[1]-T[0])/(pow(4,i+1)-1); }printf("%f\n",T[11]); }4. 用定步长四阶Runge-Kutta 求解一、程序要求用定步长四阶法求解 y1’=1y2’=y3 y3’=1000-1000y2-100y3(y1(0)=y2(0)=y3(0)=0) h=0.0005,打印yi(0.025),yi(0.045),yi(0.085),yi(0.1),(i=1,2,3)h =0.0005,打印y i (0.025) ,⎪⎪⎪⎪⎩⎪⎪⎪⎪⎨⎧===--===0)0(0)0(0)0(10010001000//1/321323321y y y y y dt dy ydt dy dt dyy i(0.045) ,y i(0.085) ,y i(0.1) ,(i=1,2,3)#include "stdafx.h"#include <stdio.h>#include <math.h>double F(double x,double y[4],double f[4]){f[1]=0*x+0*y[1]+0*y[2]+0*y[3]+1;f[2]=0*x+0*y[1]+0*y[2]+1*y[3]+0;f[3]=0*x+0*y[1]-1000*y[2]-100*y[3]+1000;return(1);}void main(){double F(double x,double y[4],double f[4]);doubleh=0.0005,x=0,Y[4],k[5][4],s[4],f[4],det,m[4]={0.025,0.045,0.085,0.1};int i,j,t; for(t=0;t<=3;t++)/*龙格-库塔算法*/{for(j=0;j<=3;j++)Y[j]=0; //每求一组值后将初值条件还原为0 for(i=1;i<=int(m[t]/h);i++){for(j=1;j<=3;j++)s[j]=Y[j];det=F(x,s,f);for(j=1;j<=3;j++)k[1][j]=h*f[j]; /*四阶古典公式中的k值和求和的计算*/ for(j=1;j<=3;j++)s[j]=Y[j]+0.5*k[1][j];det=F(x+0.5*h,s,f);for(j=1;j<=3;j++)k[2][j]=h*f[j];for(j=1;j<=3;j++)s[j]=Y[j]+0.5*k[2][j];det=F(x+0.5*h,s,f);for(j=1;j<=3;j++)k[3][j]=h*f[j];for(j=1;j<=3;j++)s[j]=Y[j]+k[3][j];det=F(x+h,s,f);for(j=1;j<=3;j++)k[4][j]=h*f[j];for(j=1;j<=3;j++)Y[j]=Y[j]+(k[1][j]+2*k[2][j]+2*k[3][j]+k[4][j])/6;x+=h;} for(j=1;j<=3;j++)printf("y[%d](%f)=%f ",j,m[t],Y[j]); printf("\n"); } } 5.⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡=40.00001 4.446782 2.213474- 0.238417 1.784317 0.037585- 1.010103- 3.123124 2.031743- 4.446782 30.719334 3.123789 1.103456- 2.121314 0.71828- 0.336993 1.112348 3.067813 2.213474- 3.123789 14.7138465 0.103458- 3.111223- 2.101023 1.056781- 0.784165- 1.7423820.238417 1.103456- 0.103458- 9.789365 0.431637 3.741856- 1.836742 1.563849 0.718719 1.7843172.1213143.111223- 0.431637 19.8979184.101011 2.031454 2.189736 0.113584-0.037585- 0.71828- 2.101023 3.741856- 4.101011 27.108437 3.123848 1.012345- 1.112336 1.010103- 0.336993 1.056781- 1.836742 2.031454 3.123848 15.567914 3.125432- 1.061074- 3.123124 1.112348 0.784165- 1.563849 2.189736 1.012345- 3.125432- 19.141823 2.115237 2.031743- 3.067813 1.742382 0.718719 0.113584- 1.112336 1.061074- 2.115237 12.38412A Tb )5.6784392- 4.719345 1.1101230 86.612343- 1.784317 0.84671695 25.173417- 33.992318 2.1874369(= 用列主元消去法求解Ax=b 。
东华大学数值分析上机样题

2003---2004 学年 第二学期 (共1页)踏实学习,弘扬正气;诚信做人,诚实考试;作弊可耻,后果自负。
课程名称 数值分析(上机部分) 使用专业研一A 卷2004年6月 时间:60分钟 每题20分班级 学号 机号 姓名考试规则:邻座试卷不同,开卷,解答全部写在考卷上;开考后不允许用软盘,也不允许上网.注:数据结果均精确到小数点后第四位;要求写出M 函数(如果需要)、MATLAB 命令和计算结果.1、在一天24小时内,从零点开始每间隔2小时测得的环境温度数据分别为 12,9,9,10,18 ,24,28,27,25,20,18,15,13, 用三次样条插值推测中午13点时的温度.(采用系统默认的边界条件) 命令:>> x=0:2:24;y=[12 9 9 10 18 24 28 27 25 20 18 15 13];>> interp1(x,y,13,'spline')答案:13点时的温度为 27.8725或>> pp=spline(x,y);>> ppval(pp,13)13点时的温度为 27.87252、利用教材程序nagauss2.m (P. 20)求解线性方程组⎪⎩⎪⎨⎧=++=++=++11163185221482321321321x x x x x x x x x .命令:>> a=[2 8 14;2 5 8 ;3 6 11];b=[2 1 1]';>> nagauss2(a,b)答案:0.33330.33330x -⎛⎫ ⎪= ⎪ ⎪⎝⎭3.求积分dx x e ⎰5)log(log 的近似值.命令:>> x=linspace(exp(1),5,100);>> y=log(log(x));>> trapz(x,y)答案:积分值为 0.63994、曲线x x y 32-=与2-=x e y 在点(0.2,-0.6)附近有交点,求其交点横坐标的近似值.命令:>> fzero('x^2-3*x-exp(x)+2',[0 2]) 答案: 横坐标为 0.25755、设30,,1,0,51 0 =+=⎰n dx x x y n n 利用关系式n y y n n 151=+-,可得如下算法:⎪⎩⎪⎨⎧=-=-=-30,,2,1,515ln 6ln 10 n y n y y n n ,试依此算法编程计算3010,,,y y y (写出302928,,y y y 即可),并说明一下此算法是否稳定. 程序:testa.mclear;format long;a=zeros(1,31);a(1)=log(6)-log(5);for n=2:31a(n)=1/(n-1)-5*a(n-1);enda命令:testa答案:28y =1841.9759 29y =-9209.8450 30y =46049.2582 算法不稳定。
数值分析上机题

数值剖析上机题习题 117.(上机题)舍入偏差与有效数N11 3 11设S N2,其精准值为 。
2 2 NN 1j 2j1(1)编制按从大到小的次序S N111,计算 S N 的通用程序。
1 32 1L22N 21(2)编制按从小到大的次序S N111,计算 S N 的通用程序。
1(N L22N 2 1)2 11(3)按两种次序分别计算 S2 ,S 4 , S 6,并指出有效位数。
(编制程序时用单精度)1010 10(4)经过本上机题你理解了什么?按从大到小的次序计算 S N 的通用程序为: 按从小到大的次序计算 S N 的通用程序为:#include<iostream.h> #include<iostream.h> float sum(float N) float sum(float N) {{float j,s,sum=0; float j,s,sum=0; for(j=2;j<=N;j++) for(j=N;j>=2;j--) {{s=1/(j*j-1); s=1/(j*j-1); sum+=s;sum+=s;}}return sum;return sum;}}从大到小的次序的值从小到大的次序的值精准值有效位数从大到小从小到大65 S 10244S 10436S 106经过本上机题, 看出按两种不一样的次序计算的结果是不同样的,按从大到小的次序计算的值与精准值有较大的偏差, 而按从小到大的次序计算的值与精准值符合。
从大到小的次序计算获得的结果的有效位数少。
计算机在进行数值计算时会出现“大数吃小数”的现象,导致计算结果的精度有所降低, 我们在计算机中进行同号数的加法时,采纳绝对值较小者先加的算法,其结果的相对偏差较小。
习题 220.(上机题) Newton 迭代法(1)给定初值x0及允许偏差,编制Newton法解方程f ( x)0 根的通用程序。
(2)给定方程 f ( x) x3 / 3 x 0 ,易知其有三个根x1 3 , x20 , x3 3 。
- 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。
如有帮助欢迎下载支持数值分析上机题姓名:陈作添学号: 040816习题 120.(上机题)舍入误差与有效数N11 3 1 1设S N,其精确值为 。
22 2 NN 1j 2j1(1)编制按从大到小的顺序111 ,计算 S 的通用程序。
S N1 321N 21 N22(2)编制按从小到大的顺序111,计算 S 的通用程序。
S N1(N 1)2 122 1NN 2(3)按两种顺序分别计算S 102 , S 104 , S 106 ,并指出有效位数。
(编制程序时用单精度)(4)通过本上机题,你明白了什么?按从大到小的顺序计算 S N 的通用程序为: 按从小到大的顺序计算 S N 的通用程序为:#include<iostream.h> #include<iostream.h> float sum(float N) float sum(float N) {{float j,s,sum=0; float j,s,sum=0; for(j=2;j<=N;j++) for(j=N;j>=2;j--) {{s=1/(j*j-1); s=1/(j*j-1); sum+=s;sum+=s;}}return sum;return sum;}}从大到小的顺序的值从小到大的顺序的值精确值有效位数从大到小从小到大0.7400490.740050.74004965 S 1020.7498520.74990.74994 4S 1040.7498520.7499990.74999936S 106通过本上机题, 看出按两种不同的顺序计算的结果是不相同的,按从大到小的顺序计算的值与精确值有较大的误差, 而按从小到大的顺序计算的值与精确值吻合。
从大到小的顺序计算得到的结果的有效位数少。
计算机在进行数值计算时会出现“大数吃小数”的现象,导致计算结果的精度有所降低,我们在计算机中进行同号数的加法时,采用绝对值较小者先加的算法,其结果的相对误差较小。
习题 220.(上机题) Newton 迭代法(1)给定初值x 0及容许误差,编制 Newton 法解方程f ( x)0 根的通用程序。
(2)给定方程f ( x) x3/ 3x0 ,易知其有三个根x 3 ,x2 0, x3 。
13 1.由 Newton 方法的局部收敛性可知存在0 ,当x0(,) 时,Newton迭代序列收敛于根 x2。
试确定尽可能大的。
2.试取若干初始值,观察当x0 (, 1),(1, ),(,) , (,1),(1, ) 时Newton 序列是否收敛以及收敛于哪一个根。
(3)通过本上机题,你明白了什么?解:( 1)编制的通用程序:#include<iostream.h>{#include<math.h>float x0,x1,a;#define eps0.000001/给定容许误差int k=0;float f(float x)//定义函数 f(x)cout<<" 请输入初值 x0:";{cin>>x0;float f;dof=x*x*x/3-x;//f(x) 的表达式 ;{return(f);a=-f(x0)/df(x0);}x1=x0+a;float df(float x)//定义函数 df(x) ,k++;计算 f(x) 的导函数x0=x1;{}float df;while(fabs(a)>eps);df=x*x-1;//f(x) 导函数的表达式 ;cout<<k<<'\t'<<x0;return (df);//输出迭代的次数和根值}}void main(void)( 2)计算迭代序列收敛于根x2的尽可能大的的函数为:#include<iostream.h>return(f);#include<math.h>}void delay(int n)// 定义延时函数float df(float x)// 定义函数 df(x) ,计算 f(x) {for(n=10000;n>0;n--);}的导函数#define eps 0.000001{float f(float x)//定义函数 f(x)float df;{df=x*x-1;//f(x) 导函数的表达式 ;float f;return (df);f=x*x*x/3-x;//f(x) 的表达式 ;}int judgement(float z)return 0;{}int count=5;void main(void)float x0,x1,type,type1;{x0=z;float delta=0;while(count-->0)int flag=1;{while(flag==1)x1=x0-f(x0)/df(x0);{type=fabs(x1);cout<<" 方程的根为 :"<<'\n';type1=fabs(x1-x0); // 调试值用delta+=eps;cout<<"count="<<count<<'\t'<<"type=flag=judgement(delta);"<<type<<'\t'<<"type1="<<type1<<'\n';}if(fabs(x1-x0)<eps)cout<<" 输出方程根收敛的区间return 1;值:\n";x0=x1;cout<<delta-eps; // 输出收敛的区间值delay(30000);//调试值用}}当步长为 0.001 时,程序计算出的δ的为δ =0.774 ,即在区间( -0.774, 0.774)内迭代序列收敛于 0。
对于不同得初始值收敛于不同的根, x0在( -∞, -1)内收敛于x1*,在( -0.774, 0.774)内收敛于x2,在( 1,+∞)内收敛于x3,但在内( 0.774,1)和(- 1,0.774)均可能收敛于x1*和 x3。
x1*, x2, x3分别为方程的精确解。
分析:对于不同的初值,迭代序列会收敛于不同的根,所以在某个区间内求根对于初值的选取有很大的关系。
产生上述结果的原因是区间不满足大范围收敛的条件。
习题 335.(上机题)列主元三角分解法对于某电路的分析,归结为求解线性方程组RI=V。
(1)编制解 n 阶线性方程组 Ax=b 的列主元三角分解法的通用程序;(2)用所编制的程序解线性方程组RI=V,并打印出解向量,保留五位有效数;(3)本编程之中,你提高了哪些编程能力?程序为:#include<iostream.h>cin>>n;#include<math.h>cout<<" 输入数组 a:"<<endl;void main(void)for(i=1;i<=n;i++){for(j=1;j<=(n+1);j++)int i,j,n,k,q;cin>>a[i][j];//给矩阵 a 赋值float a[10][11],s[10],s1[10];for(i=1;i<=n;i++)cout<<" 请输入 n 的值 :";{for(j=1;j<=(n+1);j++)cout<<a[i][j]<<'\t';cout<<'\n';}//输出数组acout<<"'''''''''''''''''''''''''"<<'\n';//进行第一行和第一列元素的求取'''''''''''''''''''''''''//int t=1;for(i=1;i<=n;i++){s[i]=a[i][1];}float max=fabs(s[1]);for(i=2;i<=n;i++)if(fabs(s[i])>max){max=fabs(s[i]);t=i;}for(j=1;j<=(n+1);j++){float b=a[1][j];a[1][j]=a[t][j];a[t][j]=b;}//进行第一列主元互换for(i=2;i<=n;i++)a[i][1]=a[i][1]/max; //第一列除以 a[1][1] for(i=1;i<=n;i++){for(j=1;j<=(n+1);j++)cout<<a[i][j]<<'\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+=a[i][q]*a[q][k];s1[i]=a[i][k]-sum;}int l=k;float m=fabs(s1[k]);for(i=k;i<=n;i++)// 比较第 k 步分解的第k 列值的大小{if(fabs(s1[i])>m){m=fabs(s1[i]);l=i;// 返回行值}}for(j=1;j<=n+1;j++)//交换两行元素{float s2=a[k][j];a[k][j]=a[l][j];a[l][j]=s2;}for(j=k;j<=n+1;j++)//算出第 k 行行元素的值{float sum1=0;for(q=1;q<k;q++)sum1+=a[k][q]*a[q][j];a[k][j]=a[k][j]-sum1;}for(i=k+1;i<=n;i++)//算出第 k 列列元素的值{float sum2=0;for(q=1;q<k;q++)sum2+=a[i][q]*a[q][k];a[i][k]=(a[i][k]-sum2)/(a[k][k]);}}//第 k 步分解结束for(i=1;i<=n;i++){for(j=1;j<=(n+1);j++)cout<<a[i][j]<<'\t';cout<<'\n';}//输出改变后的数组//输出解 '''''''''''''''''''''''''''''''''''''''''''''''''''''''//float x[10];for(i=n-1;i>=1;i--){x[n]=a[n][n+1]/a[n][n];float sum3=0;for(i=1;i<=n;i++)for(j=i+1;j<=n;j++){sum3+=a[i][j]*x[j];cout<<'x'<<i<<'='<<x[i]<<endl;x[i]=(a[i][n+1]-sum3)/a[i][i];}//输出解向量}//回代过程}结果:方程的解为:x1= -0.28923 , x2= 0.34544 ,x3= -0.71281 ,x4= -0.22061 , x5= -0.43040 , x6= 0.15431 ,x7= -0.057823 , x8= 0.20105 ,x9= 0.29023 。