0511数值分析作业一
数值分析试题集
..数值分析试题集(试卷一)一( 10 分)已知 x 1* 1.3409 ,x 2* 1.0125 都是由四舍五入产生的近似值, 判断 x 1*x 2* 及 x 1* x 2*有几位有效数字。
二( 10 分)由下表求插值多项式x 01 2 y2 34 y1- 1三( 15 分)设 f ( x)C 4 [a,b] , H ( x )是满足下列条件的三次多项式H (a) f (a) , H (b) f (b) , H (c)f (c) , H (c) f (c)( a c b )求 f (x)H ( x) ,并证明之。
12四( 15 分)计算13 dx ,10 2。
x五( 15 分)在 [0,2]上取 x 0 0 , x 1 1 , x 22 ,用二种方法构造求积公式,并给出其公式的代数精度。
六( 10 分)证明改进的尢拉法的精度是 2 阶的。
七( 10 分)对模型 yy , 0 ,讨论改进的尢拉法的稳定性。
八( 15分)求方程 x 34x 2 7x 1 0 在 -1.2 附近的近似值,10 3。
-----------------------------------------------------------------------------------------------------------------------------(试卷二)一填空( 4*2 分)1 {k ( x) } k 0 是区间 [0, 1]上的权函数为( x) x 2 的最高项系数为 1 的正交多项式族,其中10 (x)1,则x0 ( x) dx ------------------- , 1 ( x) ------------------。
2 12 A,则 A1 4----------- ,( A) ----------------- 。
a 1 2 时, A 可作 LU 分解。
3 设 A,当 a 满足条件 ---------------- 14..4 设非线性方程 f ( x) (x33x23x1)( x 3) 0 ,其根 x1* 3 , x2*1,则求 x1* 的近似值时,二阶局部收敛的牛顿迭代公式是--------------------------- 。
数值分析大作业一
数值分析大作业一一、算法设计方案1、求λ1和λ501的值:思路:采用幂法求出按模最大特征值λmax,该值必为λ1或λ501,若λmax小于0,则λmax=λ1;否则λmax=λ501。
再经过原点平移,使用幂法迭代出矩阵A-λmax I的特征值,此时求出的按模最大特征值即为λ1和λ501的另一个值。
2、求λs的值:采用反幂法求出按模最小的特征值λmin即为λs,其中的方程组采用LU分解法进行求解。
3、求与μk最接近的特征值:对矩阵A采用带原点平移的反幂法求解最小特征值,其中平移量为:μk。
4、A的条件数cond(A)=| λmax/λmin|;5、A的行列式的值:先将A进行LU分解,再求U矩阵对角元素的乘积即为A 行列式的值。
二、源程序#include<iostream>#include<iomanip>#include<math.h>#define N 501#define E 1.0e-12 //定义精度常量#define r 2#define s 2using namespace std;double a[N];double cc[5][N];void init();double mifa();double fmifa();int max(int aa,int bb);int min(int aa,int bb);int max_3(int aa,int bb,int cc);void LU();void main(){double a1,a2,d1,d501=0,ds,det=1,miu[39],lamta,cond;int i,k;init();/*************求λ1和λ501********************/a1=mifa();if(a1<0)d1=a1; //若小于0则表示λ1的值elsed501=a1; //若大于0则表示λ501的值for(i=0;i<N;i++)a[i]=a[i]-a1;a2=mifa()+a1;if(a2<0)d1=a2; //若小于0则表示λ1的值elsed501=a2; //若大于0则表示λ501的值cout<<"λ1="<<setiosflags(ios::scientific)<<setprecision(12)<<d1<<"\t";cout<<"λ501="<<setiosflags(ios::scientific)<<setprecision(12)<<d501<<endl;/**************求λs*****************/init();ds=fmifa();cout<<"λs="<<setiosflags(ios::scientific)<<setprecision(12)<<ds<<endl;/**************求与μk最接近的特征值λik**************/cout<<"与μk最接近的特征值λik:"<<endl;for(k=0;k<39;k++){miu[k]=d1+(k+1)*(d501-d1)/40;init();for(i=0;i<N;i++)a[i]=a[i]-miu[k];lamta=fmifa()+miu[k];cout<<"λi"<<k+1<<"\t\t"<<setiosflags(ios::scientific)<<setprecision(12)<<lamta<<en dl;}/**************求A的条件数**************/cout<<"矩阵A的条件式";cond=abs(max(abs(d1),abs(d501))/ds);cout<<"cond="<<setiosflags(ios::scientific)<<setprecision(12)<<cond<<endl;/**************求A的行列式**************/cout<<"矩阵A的行列式";init();LU();for(i=0;i<N;i++){det*=cc[2][i];}cout<<"det="<<setiosflags(ios::scientific)<<setprecision(12)<<det<<endl;system("pause");}/**************初始化函数,给a[N]赋值*************/void init(){int i;for(i=1;i<=501;i++)a[i-1]=(1.64-0.024*i)*sin((double)(0.2*i))-0.64*exp((double)(0.1/i)); }/**************幂法求最大绝对特征值**************/double mifa(){int i,k=0;double u[N],y[N]={0},b=0.16,c=-0.064,Beta_=0,error;for(i=0;i<501;i++)u[i]=1; //令u[N]=1for(k=1;k<2000;k++) //控制最大迭代次数为2000{/***求y(k-1)***/double sum_u=0,gh_sum_u;for(i=0;i<N;i++){sum_u+=u[i]*u[i]; }gh_sum_u=sqrt(sum_u);for(i=0;i<N;i++){y[i]=u[i]/gh_sum_u;}/****求新的uk****/u[0]=a[0]*y[0]+b*y[1]+c*y[2];u[1]=b*y[0]+a[1]*y[1]+b*y[2]+c*y[3]; //前两列和最后两列单独拿出来求中D间的循环求for(i=2;i<N-2;i++){u[i]=c*y[i-2]+b*y[i-1]+a[i]*y[i]+b*y[i+1]+c*y[i+2];}u[N-2]=c*y[N-4]+b*y[N-3]+a[N-2]*y[N-2]+b*y[N-1];u[N-1]=c*y[N-3]+b*y[N-2]+a[N-1]*y[N-1];/***求beta***/double Beta=0;for(i=0;i<N;i++){Beta+=y[i]*u[i];}//cout<<"Beta"<<k<<"="<<Beta<<"\t"; 输出每次迭代的beta /***求误差***/error=abs(Beta-Beta_)/abs(Beta);if(error<=E) //若迭代误差在精度水平内则可以停止迭代{return Beta;} //控制显示位数Beta_=Beta; //第个eta的值都要保存下来,为了与后个值进行误差计算 }if(k==2000){cout<<"error"<<endl;return 0;} //若在最大迭代次数范围内都不能满足精度要求说明不收敛}/**************反幂法求最小绝对特¬征值**************/double fmifa(){int i,k,t;double u[N],y[N]={0},yy[N]={0},b=0.16,c=-0.064,Beta_=0,error;for(i=0;i<501;i++)u[i]=1; //令u[N]=1for(k=1;k<2000;k++){double sum_u=0,gh_sum_u;for(i=0;i<N;i++){sum_u+=u[i]*u[i]; }gh_sum_u=sqrt(sum_u);for(i=0;i<N;i++){y[i]=u[i]/gh_sum_u;yy[i]=y[i]; //用重新赋值,避免求解方程组的时候改变y的值}/****LU分解法解方程组Au=y,求新的***/LU();for(i=2;i<=N;i++){double temp_b=0;for(t=max(1,i-r);t<=i-1;t++)temp_b+=cc[i-t+s][t-1]*yy[t-1];yy[i-1]=yy[i-1]-temp_b;}u[N-1]=yy[N-1]/cc[s][N-1];for(i=N-1;i>=1;i--){double temp_u=0;for(t=i+1;t<=min(i+s,N);t++)temp_u+=cc[i-t+s][t-1]*u[t-1];u[i-1]=(yy[i-1]-temp_u)/cc[s][i-1];}double Beta=0;for(i=0;i<N;i++){Beta+=y[i]*u[i];}error=abs(Beta-Beta_)/abs(Beta);if(error<=E){return (1/Beta);}Beta_=Beta;}if(k==2000){cout<<"error"<<endl;return 0;} }/**************求两数最大值的子程序**************/int max(int aa,int bb){return(aa>bb?aa:bb);}/**************求两数最小值的子程序**************/int min(int aa,int bb){return(aa<bb?aa:bb);}/**************求三数最大值的子程序**************/int max_3(int aa,int bb,int cc){ int tt;if(aa>bb)tt=aa;else tt=bb;if(tt<cc) tt=cc;return(tt);}/**************LU分解**************/void LU(){int i,j,k,t;double b=0.16,c=-0.064;/**赋值压缩后矩阵cc[5][501]**/for(i=2;i<N;i++)cc[0][i]=c;for(i=1;i<N;i++)cc[1][i]=b;for(i=0;i<N;i++)cc[2][i]=a[i];for(i=0;i<N-1;i++)cc[3][i]=b;for(i=0;i<N-2;i++)cc[4][i]=c;for(k=1;k<=N;k++){for(j=k;j<=min(k+s,N);j++){double temp=0;for(t=max_3(1,k-r,j-s);t<=k-1;t++)temp+=cc[k-t+s][t-1]*cc[t-j+s][j-1];cc[k-j+s][j-1]=cc[k-j+s][j-1]-temp;}//if(k<500){for(i=k+1;i<=min(k+r,N);i++){double temp2=0;for(t=max_3(1,i-r,k-s);t<=k-1;t++)temp2+=cc[i-t+s][t-1]*cc[t-k+s][k-1];cc[i-k+s][k-1]=(cc[i-k+s][k-1]-temp2)/cc[s][k-1];}}}}三、程序结果。
数值分析大作业
数值分析上机作业(一)一、算法的设计方案1、幂法求解λ1、λ501幂法主要用于计算矩阵的按模最大的特征值和相应的特征向量,即对于|λ1|≥|λ2|≥.....≥|λn|可以采用幂法直接求出λ1,但在本题中λ1≤λ2≤……≤λ501,我们无法判断按模最大的特征值。
但是由矩阵A的特征值条件可知|λ1|和|λ501|之间必然有一个是最大的,通过对矩阵A使用幂法迭代一定次数后得到满足精度ε=10−12的特征值λ0,然后在对矩阵A做如下的平移:B=A-λ0I由线性代数(A-PI)x=(λ-p)x可得矩阵B的特征值为:λ1-λ0、λ2-λ0…….λ501-λ0。
对B矩阵采用幂法求出B矩阵按模最大的特征值为λ∗=λ501-λ0,所以λ501=λ∗+λ0,比较λ0与λ501的大小,若λ0>λ501则λ1=λ501,λ501=λ0;若λ0<λ501,则令t=λ501,λ1=λ0,λ501=t。
求矩阵M按模最大的特征值λ的具体算法如下:任取非零向量u0∈R nηk−1=u T(k−1)∗u k−1y k−1=u k−1ηk−1u k=Ay k−1βk=y Tk−1u k(k=1,2,3……)当|βk−βk−1||βk|≤ε=10−12时,迭终终止,并且令λ1=βk2、反幂法计算λs和λik由已知条件可知λs是矩阵A 按模最小的特征值,可以应用反幂法直接求解出λs。
使用带偏移量的反幂法求解λik,其中偏移量为μk=λ1+kλ501−λ140(k=1,2,3…39),构造矩阵C=A-μk I,矩阵C的特征值为λik−μk,对矩阵C使用反幂法求得按模最小特征值λ0,则有λik=1λ0+μk。
求解矩阵M按模最小特征值的具体算法如下:任取非零向量u 0∈R n ηk−1= u T (k−1)∗u k−1y k−1=u k−1ηk−1 Au k =y k−1βk =y T k−1u k (k=1,2,3……)在反幂法中每一次迭代都要求解线性方程组Au k =y k−1,当K 足够大时,取λn =1βk 。
数值分析作业(完整版)
的逆阵 A ,用左除命令 A \ E 检验你的结果。
clc clear close all A=[1 1 1 1 1;1 2 3 4 5;1 3 6 10 15;1 4 10 20 35;1 5 15 35 70]; fprintf('对上述矩阵进行列主元素分解:\n') for i=1:1:r-1 [mx,ro]=max(abs(A(i:r,i))); % 寻找a阵第i列的最大值 [A(i,:),A(ro+i-1,:)]=exchange(A(i,:),A(ro+i-1,:)); % 进行行与行交换 for j=i+1:1:r A(j,:)=A(j,:)-A(j,i)/A(i,i)*A(i,:); end A End %--矩阵A的逆阵 A1=inv(A) %--左除验证 E=eye(5); A2=A\E % 5x5单位阵 % A阵的逆矩阵 % 输出每次交换后的A
第一章
1、计算积分 I n
Code: clc clear close all n=9; %--梯形积分法 x=0:0.01:1; y=(x.^n).*exp(x-1); In = trapz(x,y); In2=vpa(In,6) % 6位有效数字 %--高精度积分法 F = @(x1)(x1.^n).*exp(x1-1); s = quad(F,0,1); s1=vpa(s,6)
0
0, 0, 0, 0, 0 。
T
if abs(er(:,i-1))<=e fprintf('在迭代 %d 次之后,满足精度要求,x向量的值如下:\n',i); fprintf('x1=%.5f, x2=%.5f, x3=%.5f, x4=%.5f, x5=%.5f\n',x(1,i),x(2,i),x(3,i),x(4,i),x(5,i)); break end end %--绘图 figure(1) plot(1:1:i,x(1,:),'b',1:1:i,x(2,:),'k',1:1:i,x(3,:),'g',1:1:i,x(4,:), 'r',1:1:i,x(5,:),'c') legend('x1','x2','x3','x4','x5') grid on title('Jacobi迭代法——x值随迭代次数变化曲线') figure(2) plot(1:1:i-1,er(1,:),'b',1:1:i-1,er(2,:),'k',1:1:i-1,er(3,:),'g',1:1: i-1,er(4,:),'r',1:1:i-1,er(5,:),'c') legend('△x1','△x2','△x3','△x4','△x5') grid on title('Jacobi迭代法——△x值随迭代次数变化曲线') %% fprintf('\n-------------Gauss-Seidel迭代法---------------------\n'); U=-(A1-D); L=-(A2-D); DL_1=inv(D-L); M1=DL_1*U; b2=DL_1*b; x1(:,1)=M1*x0+b2; for j=2:1:100 x1(:,j)=M1*x1(:,j-1)+b2; er1(:,j-1)=x1(:,j)-x1(:,j-1); if abs(er1(:,j-1))<=e fprintf('在迭代 %d 次之后,满足精度要求,x向量的值如下:\n',j); fprintf('x1=%.5f, x2=%.5f, x3=%.5f, x4=%.5f, x5=%.5f\n',x1(1,j),x1(2,j),x1(3,j),x1(4,j),x1(5,j)); break end end %--绘图 figure(3) plot(1:1:j,x1(1,:),'b',1:1:j,x1(2,:),'k',1:1:j,x1(3,:),'g',1:1:j,x1(4 ,:),'r',1:1:j,x1(5,:),'c') legend('x1','x2','x3','x4','x5')
(完整版)数值分析第一次作业
问题1:20.给定数据如下表:试求三次样条插值S(x),并满足条件 (1)S`(0.25)=1.0000,S`(0.53)=0.6868; (2)S ’’(0.25)=S ’’(0.53)=0。
分析:本问题是已知五个点,由这五个点求一三次样条插值函数。
边界条件有两种,(1)是已知一阶倒数,(2)是已知自然边界条件。
对于第一种边界(已知边界的一阶倒数值),可写出下面的矩阵方程。
⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡432104321034322110d M M M M M 200020000020022d d d d λμμλμλμλ其中μj =j1-j 1-j h h h +,λi=j1-j j h h h +,dj=6f[x j-1,x j ,x j+1], μn =1,λ0=1对于第一种边界条件d 0=0h 6(f[x 0,x 1]-f 0`),d n =1-n h 6(f`n-f `[x n-1,x n ]) 解:由matlab 计算得:由此得矩阵形式的线性方程组为:⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎣⎡ 2.1150-2.4286-3.2667-4.3143-5.5200-M M M M M 25714.00001204286.000004000.026000.0006429.023571.0001243210解得 M 0=-2.0286;M 1=-1.4627;M 2= -1.0333; M 3= -0.8058; M 4=-0.6546S(x)=⎪⎪⎩⎪⎪⎨⎧∈-+-+-∈-+-+-∈-+-+-∈-+-+-]53.0,45.0[x 5.40x 9.1087x 35.03956.8.450-x 1.3637-x .5301.67881- ]45.0,39.0[x 9.30x 11.188x 54.010.418793.0-x 2.2384-x .450(2.87040-]39.0,30.0[x 03.0x 6.9544x 9.30 6.107503.0-x 1.9136-x .3902.708779-]30.0,25.0[x 5.20x 10.9662x 0.3010.01695.20-x 4.8758-x .3006.76209-33333333),()()()(),()()()),()()()(),()()()(Matlab 程序代码如下:function tgsanci(n,s,t) %n代表元素数,s,t代表端点的一阶导。
2011数值分析试题及答案
解:由 x( k 1) Mx( k ) g 和 x* Mx* g 可得:
x( k 1) x* M ( x( k ) x* ) , k 0,1,2,...
递推的: x( k ) x* M k ( x(0) x* ) 设 y 是矩阵 M 属于特征值 的特征向量,取 x(0) y x* ,则有:
1 1
1 0 0 0 0.3 0.2 0 0.3 0.2 0 0. 4 G ( D L) U 0 1 0 0 0 0.4 0 1 0 1 0 0 0 0 0. 3 0. 2
3 5x f ( xk ) xk a a , xk 1 k 2 , k 0,1,2,... xk xk 2 6 f ( xk ) 6 xk 6 xk
一、解答下列各题: (每题 5 分,共 30 分) 1.设近似值 x 具有 5 位有效数字,则 x 的相对误差限为多少? 解:记 x 0.a1a2 ...10 ,则 x 的相对误差为:
五、 (4 分)设矩阵 M 是 n 阶方阵, M 有一个绝对值小于 1 的特征值 ,且方程 组 x Mx g 有 唯 一 解 x * , 证 明 : 存 在 初 始 向 量 x ( 0 ) 使 迭 代 格 式 :
x ( k 1) Mx ( k ) g , k 0,1,2,...产生的序列 {x ( k ) } 收敛到 x * .
数值分析作业
数值分析课后作业:习题一1.在字长为3的十进制计算机上计算f (3.33)和g (3.33),其中f(x)=x 4-x 3+3x 2+x-2,g(x)=(((x-1)x+3)x+1)x-2解: m=3; f=@(x)digit(digit(x^4,m)- digit(x^3,m)+ digit(3*x^2,m)+ digit(x-2,m),m); g=@(x)digit(digit(digit( digit(digit(digit( (x-1)*x,m)+3,m)*x,m)+1,m)*x,m)-2,m); f(3.33) g(3.33) 有ans = 121 ans =121 2.下列各近似值的绝对误差限都是1021⨯-3,试指出它们各有几位有效数字:x=1.00052, y=0.05, z=0.00052.解:当 x=1.00052时, 由丨X*—X 丨 ≤0.5×10-3 得 x=1.00052 有四位有效数字; 同理 y=0052 有两位有效数字 Z=0.00052有零位有效数字 3,计算圆的面积,要使其相对误差限为1%,问测量半径r 允许的相对误差限是多少? 解:设圆的面积为S , 由题意有|e(S)|≤1%。
又S=πr 2 dS=2πr dr 所以 dS/S=(2πrdr)/(πr 2)=2(dr/r)∴|e(r)|≈21|e(S)|≤0.5×1%=0.5% 11.数组与矩阵是Matlab 编程的基础,试学习Matlab 的数组与矩阵的表示方法,并举例介绍数组、矩阵的常见运算. 解:>> syms a b c d; >> a=[1 2 3];>> b=[4 5 6];>> a+bans =5 7 9>> b-aans =3 3 3>> a.*bans =4 10 18 >> a.^2 ans = 1 4 9>> c=[1 2 3;1 2 3;1 2 3];>> d=[4 5 6;4 5 6;4 5 6];>> cc = 1 2 3 1 2 3 1 2 3d = 4 5 6 4 5 6 4 5 6 >> c+dans =5 7 9 5 7 9 5 7 9>> d-cans = 3 3 33 3 33 3 3 12.学习使用Matlab 命令help 和doc 学习自己感兴趣的Matlab 的运算、函数或命令的用法,并对于任意给定的实数a,b,c,编写Matlab 程序求方程ax 2+bx+c=0的根. 解:x 1=a ac b b b 24)sgn(2---, x 2=1ax c1 x>0 其中 sgn = 0 x=0 -1 x<0 disp('Please input the coefficients of');disp('quadratic equation ax^2+bx+c=0, respectively') a=input('a='); b=input('b='); c=input('c=');m=3; if abs(a)<eps & abs(b)<eps error End if abs(a)<eps disp('Since a=0, quadrtic equation degen erates into a linear equation.') disp('The only solution of the linear equtio n is')x=digit(-c/b,m) return Enddelta=b^2-4*a*c; temp=sqrt(delta); x 1=(-b+temp)/(2*a) ; x 2=(-b-temp)/(2*a) ;err1=abs(a*x 1^2+b*x 1+c) ; err2=abs(a*x 2^2+b*x 2+c) ; if b>0x 1=(-b-temp)/(2*a) End if b<0x 1=(-b+temp)/(2*a) End if b=0x 1=temp/(2*a) Endx 2=c/(a*x 1)err1=abs(a*x 1^2+b*x 1+c) err2=abs(a*x 2^2+b*x 2+c) if abs(a)<epsdisp('Since a=0, quadrtic equation degen erates into a linear equation.')disp('The only solution of the linear equtio n is')x=digit(-c/b,m) return Enddelta=digit(digit(b^2,m)-digit(4*digit(a*c,m),m),m);temp=digit(sqrt(delta),m);x 1=digit(digit(-b+temp,m)/digit(2*a,m),m); x 2=digit(digit(-b-temp,m)/digit(2*a,m),m); err1=abs(a*x 1^2+b*x 1+c); err2=abs(a*x 2^2+b*x 2+c); if b>0x 1=digit(digit(-b-temp,m)/digit(2*a,m),m) ; End if b<0x 1=digit(digit(-b+temp,m)/digit(2*a,m),m); End if b=0x 1=digit(temp/digit(2*a,m),m); Endx 2=digit(digit(c/a,m)/x1,m) ; err1=abs(a*x 1^2+b*x 1+c) ; err2=abs(a*x 2^2+b*x 2+c) ; 14分别利用ln (1+x)=11,)1(11≤<--+∞=∑x nx nn n 和ln11...),12...53(2111253<<-++++++=-++x n x x x x x x n ,给出计算ln2的近似方法,编写相应的Matlab 程序,并比较算法运行情况. 解:方法一: x=1; s=0;for k=1:100s=s+(-1)^(k+1)*(x^k)/k; end sq=log(2)err=abs(t-q) ans= t =0.6882 q =0.6931 err = 0.0050方法二x=1/3; s=0;for k=1:2:100 s=s+(x^k)/k; end t=2*s q=log(2)err=abs(t-q) Ans= t =0.6931 q =0.6931 err =2.2204e-16所以方法二较方法一好。
数值分析习题(含标准答案)
]第一章 绪论姓名 学号 班级习题主要考察点:有效数字的计算、计算方法的比较选择、误差和误差限的计算。
1若误差限为5105.0-⨯,那么近似数有几位有效数字(有效数字的计算) 解:2*103400.0-⨯=x ,325*10211021---⨯=⨯≤-x x 故具有3位有效数字。
2 14159.3=π具有4位有效数字的近似值是多少(有效数字的计算) 解:10314159.0⨯= π,欲使其近似值*π具有4位有效数字,必需!41*1021-⨯≤-ππ,3*310211021--⨯+≤≤⨯-πππ,即14209.314109.3*≤≤π即取( , )之间的任意数,都具有4位有效数字。
3已知2031.1=a ,978.0=b 是经过四舍五入后得到的近似值,问b a +,b a ⨯有几位有效数字(有效数字的计算)解:3*1021-⨯≤-aa ,2*1021-⨯≤-b b ,而1811.2=+b a ,1766.1=⨯b a 2123****102110211021)()(---⨯≤⨯+⨯≤-+-≤+-+b b a a b a b a故b a +至少具有2位有效数字。
2123*****10210065.01022031.1102978.0)()(---⨯≤=⨯+⨯≤-+-≤-b b a a a b b a ab 故b a ⨯至少具有2位有效数字。
4设0>x ,x 的相对误差为δ,求x ln 的误差和相对误差(误差的计算)~解:已知δ=-**xx x ,则误差为 δ=-=-***ln ln xx x x x则相对误差为******ln ln 1ln ln ln xxx x xxx x δ=-=-5测得某圆柱体高度h 的值为cm h 20*=,底面半径r 的值为cm r 5*=,已知cm h h 2.0||*≤-,cm r r 1.0||*≤-,求圆柱体体积h r v2π=的绝对误差限与相对误差限。
(误差限的计算)解:*2******2),(),(h h r r r h r r h v r h v -+-≤-ππ绝对误差限为πππ252.051.02052)5,20(),(2=⨯⋅+⨯⋅⋅⋅≤-v r h v相对误差限为%420120525)5,20()5,20(),(2==⋅⋅≤-ππv v r h v 6设x 的相对误差为%a ,求nx y =的相对误差。
数值分析作业答案
数值分析作业答案【篇一:《数值分析》作业参考答案 2】>一. 选择题1. a; 2.b; 3.b; 4.d; 5.c; 6.d; 7.c; 8.b; 9.d; 10.c; 11.b;12.a; 13.a; 14.c; 15.a; 16.b; 17.d; 18.a19.d,20.c,21.a,22.d,23.c,24.c. 二. 填空题1. 3,3,3 ;2. 1,2/3;3. 100!2^100 ;4. (-1 ,1);5.x?cosx1; 6. g(x)?x?,2;1?sinx(x0?x1)(x0?x2)?(x0?xn)7. (-1, 1); 8. x; 9. 4 ; 10. 5,9 ; 11. xn?1?31c(xn?), 2; 2xnb?a1x2?112. x?x?1 ; 13. 10/9, 4; 14. 10, 55, 550; 15. b?16.17.33318.?三.1. p(x)?x2?2x?12. p(x)?f(x)3. 2.a?c?a?b11?i1?3i219.2x?x20. a? . ,,324429121016, 代数精确度为5 ,b?,a?9953.证明:cond(aa)?||aa||?||(aa)?1|| 设??max{aa的特征值的模},??max{(a?1)a?1的特征值的模},则上式=????||a||2?||a?1||2?(cond(a))2?1???24.1. (12分)l???????2?1??1???????1????0?2u?, , x??????4?1?1???????0??1??????1?31?2. (8分)seidel收敛,因为a 实正定对称阵. 迭代格式(k)?x1(k?1)?(2?x2)/2?(k?1)(k)?(x1(k?1)?x3)/2?x2?(k?1)(k?1)(k)x?(?2?x?x24)/2?3(k?1)(k?1)?x4?(1?x3)/2?12??35. p(x)??2x?1,余项cosx?p2(x)?|x(x?|?6.?6254*6426. 证明:当a?0时,结论显然成立;|ytax|当a?0时,因|yax|?||y||2||a||2||x||2,故 sup?||a||2;x?0,y?0||x||2||y||2t又aa是实对称矩阵,故存在正交阵p?(p1,p2,?,pn)t??1???tt?使得paap?d???, pi是特征值?i对应的特征向量。
数值分析报告作业第一次
第二章 习题20.试求三次样条插值S(x),并满足条件 (1)S`(0.25)=1.0000,S`(0.53)=0.6868;分析:已知两端的一阶导数值为第一种边界条件。
可写成矩阵:⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡432104321034322110d M M M M M 200020000202002d d d d λμμλμλμλ 其中μj =j1-j 1-j h h h +,λi=j1-j j h h h +,dj=6f[x j-1,x j ,x j+1], μn =1,λ0=1对于第一种边界条件d 0=0h 6(f[x 0,x 1]-f 0`),d n =1-n h 6(f`n-f `[x n-1,x n ] 解:由matlab 计算得: 、由此得矩阵形式的线性方程组为:⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎣⎡-----=⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎣⎡1150.24286.22667.33143.45200.5M M M M M 0000.25714.000010000.204286.000004000.00000.26000.0006429.00000.23571.000000.10000.243210解得M 0=-2.0286 ;M 1=-1.4627;M 2= -1.0333; M 3= -0.8058; M 4=-0.6543S(x)=⎪⎪⎩⎪⎪⎨⎧∈-+-+-∈-+-+--∈-+-+--∈-+-+--]53.0,45.0[x 5.40x 6708.0x 35.07452.0.450-x 0.4029-x .5300.3151]45.0,39.0[x 9.30x 62451.0x 54.00.800493.0-x .51670x .4506320.0]39.0,30.0[x 03.0x 5477.0x 9.309127.003.0-x 7314.0x .3907952.0]30.0,25.0[x 5.20x 5000.0x 0.300000.15.20-x 0143.1x .3008863.133333333),()()()(),()()()(),()()()(),()()()( matlab 源程序>>x=[0.25 0.30 0.39 0.45 0.53];y=[0.5 0.5477 0.6245 0.6708 0.7280];S0=1;S5=0.6868;for j=1:1:4h(j)=x(j+1)-x(j);endfor j=2:1:4r(j)=h(j)/(h(j-1)+h(j));endr(1)=1;for j=1:1:3u(j)=h(j)/(h(j)+h(j+1));endu(4)=1;for j=1:1:4f(j)=(y(j+1)-y(j))/h(j);endd(1)=6*(f(1)-S0)/h(1);d(5)=6*(S5-f(4))/h(4);for j=2:1:4d(j)=6*(f(j)-f(j-1))/(h(j-1)+h(j));enda=zeros(5,5);for i=1:1:5a(i,i)=2;endfor i=1:1:4a(i+1,i)=u(i);a(i,i+1)=r(i);endb=inv(a);M=b*d';s=csape(x,y,'complete',[1 0.6868])>> fnplt(s,'r')>> xlabel('x')>> ylabel('y')title('三次样条插值函数')plot(x,y,'o',x,y,'')>> s.coefs三次样条插值函数图xy(2)S``(0.25)=S``(0.53)=0.分析:已知两端的二阶导数只为零,可以利用自然边界条件。
数值分析作业1
拉格朗日插值法摘要:本篇综述是从插值法的原理入手,通过线性插值(一次插值) 、抛物线插值(二次插值)的分析,从特殊到一般,从简单到复杂,引入拉格朗日插值多项式。
关键词:插值法线性插值抛物插值拉格朗日插值一、引言在科学研究和其他许多实际问题中,常常有函数不便于处理和计算的情形。
有时候函数关系没有明显的解析表达式,需要根据实验观测或其他方法来确定与自变量的某些值相对应的函数值;有时函数虽有解析表达式,但是使用很不方便。
因此,希望对这些问题中的函数建立一个简单的便于计算和处理的近似表达式,即用一个简单的函数来近似代替这些不变处理的函数。
与用近似数代替准确数一样,这也是数值计算方法中最基本的概念和方法之一,也就是插值法。
二、插值法的基本原理1. 插值法设函数y =f(x)定义在区间〔a , b ]上,x°,X i,…,X n是la , b 1上取定的n+1个互异的节点,且在这些点处的函数值f (X o),f (xj,…,f (X n)为已知,即y i二f(X i),若存在一个f(x)的近似函数P(X),满足p(Xi)= yi i = 0, 1,2, , , n。
( 1),则称p( x )为f (x)的插值函数。
(1)式称为插值条件,f ( x )称为被插函数,[a , b ]称为插值区间,求p ( x )的方法就是插值法。
插值函数p ( x )在n+1个互异插值节点X i ( i = 0, 1, 2, , , n)处与f (x )相等,在其他点x就用p( x )的值作为f (x)的近似值。
这一过程称为插值,点x称为插值点。
换句话说,插值就是根据被插函数给出的函数表“插出”所要点的函数值。
用p ( X )的值作为f (x)的近似值。
不仅希望p ( X )能较好的逼近f(X),而且还希望它简单。
由于代数多项式计算既简单又便于计算,这是我们选择用多项式作为插值函数的原因。
我们数值分析课主要学习了代数多项式插值。
数值分析第一次作业答案
作业1.用如下数值表构造不超过3次的插值多项式2. P55 11题.给出概率积分⎰-=xxdxey 022π的数据表用2次插值计算,试问:(1) 当x = 0.472时,积分值等于多少? (2) 当x 为何值时,积分值等于0.5? 解:(1) 取x 0 = 0.47, x 1 = 0.48, x 2 = 0.4980.4955530040.04093346-80.1809899240.355496540.51166830.50274980.4937452=+=----⨯+----⨯+----⨯==----+----+----≈)48.049.0)(47.049.0()48.0472.0)(47.0472.0()49.048.0)(47.048.0()49.0472.0)(47.0472.0()49.0472.0)(48.047.0()49.0472.0)(48.0472.0()472.0())(())(())(())(())(())(()472.0(2120210221120121210y Lxx xx xxy x x x x x x y x x x x x x x x x x x x y(2)90.4769359350.05272367-80.4362204360.093439170.50274980.51166830.49374520.51166830.50274980.49374520.49 0.51166830.50274980.49374520.50274980.51166830.49374520.48 0.51166830.49374520.50274980.49374520.51166830.50274980.47=+=----⨯+----⨯+----⨯==----+----+----≈))(()5.0)(5.0())(()5.0)(5.0())(()5.0)(5.0()5.0())(())(())(())(())(())(()5.0(212210221120121210Lyyy y yy xy y yy y y xyyyy yy xy y y y y y x3. 证明方程e x +10x -2=0在区间[0,1]内有一个根,如果使用二分法求该区间内的根,且误差不超过10-6,试问需要二分区间[0,1]多少次?4. 设x t =451.01为准确值,x a =451.023为x t 的近似值,试求出x a 有效数字的位数及相对误差 作业答案1.解:N 2(x ) = f (0)+f [0,1](x -0)+ f [0,1,2](x -0) (x -1) 1+1×(x -0) +3×(x -0) (x -1)=3x 2-2x +1 为求得P 3(x ),根据插值条件知,P 3(x )应具有下面的形式 P 3(x )=N 2(x )+k (x -0) (x -1) (x -2),这样的P 3(x )自然满足:P 3(x i )= f (x i )由P 3’(1 )=3P 3’(1 )= N 2’(1 )+k (1-0) (1-2) =N 2’(1 )-k = 4-k=3∴ k =1∴ P 3(x )=N 2(x )+ (x -0) (x -1) (x -2)=x 3+1 3. 证明 令f (x )=e x +10x -2,∵ f (0)=-1<0,f (1)=e+8> 0∴ f (x )= e x +10x -2 =0在[0,1]有根。
数值分析第1次上机作业
《数值分析》第1次上机作业姓名:学号:2015.11.051 算法设计思路和方案因为题目只要求求解部分特征值,及最大的特征值和最小的特征值以及距离某些给定实数最近的特征值,而矩阵A 是实对称矩阵,所有特征值均为实数,故我们可以用幂法和反幂法结合适当的平移量求解。
1.1 算法分析1.1.1 幂法求解1λ,501λ我们知道幂法适用于求解矩阵 A 的绝对值最大的特征值,设矩阵A 的特征值为12n λλλ≥≥≥。
对于情况12λλ>,幂法是适用的,并可以求出特征值1λ。
然而对于12=λλ时,幂法只对12λλ=适用(即当绝对值最大的特征值是重根的情况,幂法是适用的);对于12λλ=-,绝对值最大的特征值为互为相反数的情况并不适用。
所以在用幂法求解1λ、501λ时,我们分两种情况讨论。
(1)1501λλ≠-此时,由于12501λλλ≤≤≤,我们知道矩阵A 绝对值最大的特征值必然为1λ、501λ中的一个,而且幂法适用,通过幂法迭代一定次数后就可以得到满足精度要求的特征值λ。
这时对A 做平移B A I λ=-,然后对矩阵B 用幂法。
由线性代数可知,矩阵B 的特征值为1λλ-,2λλ-,⋯,501λλ-。
若1λλ=,则有125010λλλλλλ=-≤-≤≤-,对B 用幂法即可求出其绝对值最大的特征值501'λλλ=-,随即可以得到A 的特征值501'λλλ=+。
对于另外一种情况501λλ=,则有125010λλλλλλ-≤-≤≤-=,对B 用幂法即可求出其绝对值最大的特征值1'λλλ=-,同样随即可以得到A 的特征值1'λλλ=+。
综上,对于1501λλ≠-,只需要对A 用幂法即可得到A 的绝对值最大的特征值λ,然后做平移B A I λ=-,对B 用幂法即可求出其绝对值最大的特征值'λ,那么有{}{}1501,',λλλλλ+=。
(2)1501λλ=-这种情况不能直接对矩阵A 使用幂法。
2011数值分析作业
第一章作业求积分function qiujifen = qiuzhi()y = 1.^9.*exp(1-1);y1 = 0.5*y;y2 = y1./2+1./2.*0.5.^9.*exp(0.5-1);S1 = 4.*y2./3-y1./3;for k = 2:infz = 0;for i = 2^kfor j = 1:2:i-1z = z+j./i.^9.*exp(j./i-1);endy1 = y2;y2 = y1./2+1./i.*z;S2 = 4.*y2./3-y1./3;endif max(abs(S2-S1))<15*0.5e-6 %zhong zhi tiaojian breakelseS1 = S2;endenddisp('结果为')y = y2Mesh命令画出二元函数function mesh = mesh(x,y,z)x = linspace(-10,10,100);y = x;[a,b] = meshgrid(x, y);z = exp(-abs(a))+cos(a+b)+1./(a.^2+b.^2+1); surf(x,y,z);第二章作业列主元三角分解,求逆并验证function pascal = pascal()A = [1 1 1 1 1;1 2 3 4 5;1 3 6 10 15;1 4 10 20 35;1 5 15 35 70]; n = length(A);P = eye(n);L = eye(n);U = zeros(n,n);for j = 1:nmax = j;for i = j+1:nB = zeros(1,n-j+1);%%%%%%%求最大列主元所在行数%%%%%%if j == 1 & abs(A(i,j))>abs(A(max,j))max = i;elseif j ~= 1w = abs(A(i,j)-A(i,1:j-1)*A(1:j-1,j));v = abs(A(max,j)-A(max,1:j-1)*A(1:j-1,j))if w>vmax = i;endend%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%endB = P(j,:);P(j,:) = P(max,:);P (max,:) = B;temp = A(j,j:n);A(j,j:n) = A(max,j:n);A(max,j:n) = temp;if j == 1A(j,j) = A(j,j)elseif j ~= 1A(j,j) = A(j,j)-A(j,1:j-1)*A(1:j-1,j);endendU(j,j) = A(j,j);for i = j+1:nif j == 1A(i,j) = A(i,j)./A(j,j);elseA(j,i) = A(j,i)-A(j,1:j-1)*A(1:j-1,i);A(i,j) = (A(i,j)-A(i,1:j-1)*A(1:j-1,j))./A(j,j);endU(j,i) = A(j,i);L(i,j) = A(i,j);endenddisp('列主元三角分解为:')P = PL = LU = UV = U\eye(n);V = V*(L\eye(n));V = V*P;disp('矩阵A的逆矩阵:')V = VA = [1 1 1 1 1;1 2 3 4 5;1 3 6 10 15;1 4 10 20 35;1 5 15 35 70]; disp('左除命令验证:')V = A\eye(n)运行结果:列主元三角分解为:可以验证列主元三角分解法正确求逆矩阵函数function ni = nijuzhen()A = [1 1 1 1 1;1 2 3 4 5;1 3 6 10 15;1 4 10 20 35;1 5 15 35 70]; n = length(A);U = A;E = eye(n);V = E;for i=1:nL(i,i) = 1;for j = i+1:nL(j,i) = U(j,i)/U(i,i);U(j,:) = U(j,:)-L(j,i)*U(i,:);V(j,:) = V(j,:)-L(j,i)*V(i,:);endendfor j = 1:nfor i = j:n-1,m = U(n-i,n-j+1)/U(n-j+1,n-j+1);U(n-i,:) = U(n-i,:)-m*U(n-j+1,:);V(n-i,:) = V(n-i,:)-m*V(n-j+1,:);endm = U(n-j+1,n-j+1);U(n-j+1,:) = U(n-j+1,:)/m;V(n-j+1,:) = V(n-j+1,:)/m;enddisp('矩阵A的逆矩阵:')Vdisp('左除法验证:')V = A\E第三章作业Jacobi迭代法function Jacobi = Jacobi()A = [10 1 2 3 4;1 9 -1 2 -3;2 -1 7 3 -5;3 2 3 12 -1;4 -3 -5 -1 15];b = [12 -27 14 -17 12];x = [0 0 0 0 0];n = length(x);k = 1;while 1for i = 1:ny(i) = b(i);for j = 1:nif j ~= iy(i) = y(i)-A(i,j)*x(j);endendy(i) = y(i)/A(i,i);endif norm(abs(y-x))<1e-5break;endx = y;k = k+1;enddisp('Jacobi迭代结果为:')yJacobi迭代结果为:y = 1.0000 -2.0000 3.0000 -2.0000 1.0000Gauss-Seidel迭代法function Gauss_Seidel = G_S()A = [10 1 2 3 4;1 9 -1 2 -3;2 -1 7 3 -5;3 2 3 12 -1;4 -3 -5 -1 15];b = [12 -27 14 -17 12];x = [0 0 0 0 0];n = length(x);while 1y(1) = b(1);y(n) = b(n);for j = 2:ny(1) = y(1)-A(1,j)*x(j);endy(1) = y(1)/A(1,1);for i = 2:n-1y(i) = b(i);for j = 1:i-1y(i) = y(i)-A(i,j)*y(j);endfor j = i+1:ny(i) = y(i)-A(i,j)*x(j);endy(i) = y(i)/A(i,i);endfor j = 1:n-1y(n) = y(n)-A(n,j)*y(j);endy(n) = y(n)/A(n,n);if norm(abs(y-x))<1e-5break;endx = y;enddisp('Gauss-Seidel迭代结果为:')yGauss-Seidel迭代结果为:y = 1.0000 -2.0000 3.0000 -2.0000 1.0000function tezhengzhi = tezhengzhi()A = [12 6 -6;6 16 2;-6 2 16];x = [1;1;1];%%%%%%幂法迭代三次%%%%%%for i = 1:3y = x/max(abs(x));[i j] = max(abs(x));y = x/x(j);x = A*y;end%%%%%%%%%%%%%%%%%%%%%%[i j] = max(abs(x)); %寻找绝对值对大分量的位置以确定其符号y = x/x(j);l = fix(x(j));A = A-l*eye(3);u = l;while(1)x = inv(A)*y;[i j] = max(abs(x));y = x/x(j);if abs(1/l-1/x(j))<1e-10breakelse[i j] = max(abs(x));l = x(j);endend[i j] = max(abs(x));disp('近似值为:')u+1/x(j)近似值为:18求插值多项式,计算误差,画图function newton = newton()x = linspace(-5,5,11);y = 1./(1+4.*x.^2);%%%%%%%%%%计算差商%%%%%%%%%%%%%y(i)是(i-1)阶差商%%%%%for i = 2:11for j = 11:-1:iy(j) = (y(j)-y(j-1))/(x(j)-x(j-i+1));endend %%%%%%%%%%%%%%%%%%%%%%%%%syms XK(1) = X-x(1);for i = 2:length(x)K(i)=(X-x(i)).*K(i-1);endN = y(1);for i=2:length(x)N = N+K(i-1).*y(i);enddisp('牛顿插值多项式为:')N = simple(simple(N)) %写出插值多项式%%计算99个节点近似值z = -5:0.1:5;Y = 1./(1+4.*z.^2);figure(1)plot(z,Y,'r')hold onu = y(11);m = length(z);for j=1:mfor i=10:-1:1u=y(i)+u*(z(j)-x(i));v(j)=u;endu=y(11);enddisp('节点近似值为:')vplot(z,v,'b')disp('各节点误差为:')Y-vfigure(2)plot(z,Y-v)原函数和插值多项式图形误差曲线用最小二乘法作出拟合曲线图function fit = fitting()x = [2 3 5 6 7 9 10 11 12 14 16 17 19 20];y = [106.42 108.26 109.58 109.50 109.86 110.00 109.93 110.59 110.60 110.72 110.90 110.76 111.10 111.30];xx = 1./x;yy = 1./y;c = [length(xx) sum(xx);sum(xx) sum(xx.*xx)];d = [sum(yy);sum(xx.*yy)];b = (d(2)-d(1)*c(2)/c(1))/(c(4)-c(3)*c(2)/c(1))a = (d(1)-b*c(3))/c(1)plot(x,y,'*'); % 在图上画出原函数真实点hold onx = 2:0.1:20;y = 1./(a+b./x);plot(x,y,’r’); %作出拟合曲线图考钮螺线function kaoniu = kaoniu()%%%%求x(s)关于t的积分%%%%s = -5:0.1:5;t = cos(0)+cos(0.5.*s.*s);%%%%使用simpson公式的逐次半分算法%%%%x1 = 0.5*t;x2 = x1./2+s./2.*cos(0.5.*0.5.*0.5.*s.*s);SX1 = 4.*x2./3-x1./3;z = 0;for k = 2:inffor i = 2^kfor j = 1:2:i-1z = z+cos(0.5.*s.*s.*j.*j./i./i);endx1 = x2;x2 = x1./2+s./i.*z;SX2 = 4.*x2./3-x1./3;endif max(abs(SX2-SX1))<15e-6 %终止条件breakelseSX1 = SX2;endend %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%求y(s)关于t的积分%%%%t = sin(0)+sin(0.5.*s.*s);y1 = 0.5*t;y2 = y1./2+s./2.*sin(0.5.*0.5.*0.5.*s.*s);SY1 = 4.*y2./3-y1./3;z = 0;for k = 2:inffor i = 2^kfor j = 1:2:i-1z = z+sin(0.5.*s.*s.*j.*j./i./i);endy1 = y2;y2 = y1./2+s./i.*z;SY2 = 4.*y2./3-y1./3;endif max(abs(SY2-SY1))<15e-6 %zhong zhi tiaojian breakelseSY1 = SY2;endendplot(SX2,SY2);第八章作业解非线性方程function solution = solution_to_nonline()b = 513/0.6651;a = -b;% a<x<bdelta = 0.5e-9;x = a+delta;i = 1;y = log((513+0.6651*x)/(513-0.6651*x))-x/1400/0.0198; %%%%%%%%从左断点开始判断,选取第一个非解点%%%%%%%%%%%while(1)if x >= bbreakendif y == 0s(i) = x;i = i+1;x = x+0.1;y = log((513+0.6651*x)/(513-0.6651*x))-x/1400/0.0198;%y是所取小区间的左端点函数值elsex = x+0.1;breakendend %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%while(1)%%%%%%%%%%%%%%%%%%寻找有根区间%%%%%%%%%%%%%%%%%%%while(1)y1 = log((513+0.6651*x)/(513-0.6651*x))-x/1400/0.0918;%y1是所取小区间右端点函数值if y*y1<0breakelseif y*y1 == 0s(i) = x;i = i+1;x = x+0.1;y = log((513+0.6651*x)/(513-0.6651*x))-x/1400/0.0918; x = x+0.1;elsex = x+0.1;endend%%%%%%%%%%%%%%%%%%%%%对分区间法%%%%%%%%%%%%%%%%%%if x>=513/0.6651&length(s) ~= 0x = sbreakelseif x>=513/0.6651&length(s) == 0disp('原方程无解!')endm = x-0.1;n = x;while(1)if (n-m)/2<deltas(i) = (m+n)/2;i = i+1;y = log((513+0.6651*x)/(513-0.6651*x))-x/1400/0.0918;x = x+0.1;breakelsex0 = (m+n)/2;y0 = log((513+0.6651*x0)/(513-0.6651*x0))-x0/1400/0.0198;if y1*y0 == 0s(i) = x0;i = i+1;breakelseif y1*y0>0n = x0;elsem = x0;endendendendx = -767.4126 0.0000 767.3874第九章作业RK法function RK4 = RK4()x = 0:pi/8:pi;h = pi/8;y = x;y(1) = 1;% f(x,y)是y的导数函数for i = 1:length(x)-1K1 = f(x(i),y(i));K2 = f(x(i)+0.5*h,y(i)+0.5*K1*h);K3 = f(x(i)+0.5*h,y(i)+0.5*K2*h);K4 = f(x(i)+h,y(i)+K3*h);y(i+1) = y(i)+h*(K1+2*K2+2*K3+K4)/6;endy %输出RK算法解y1 = cos(x)+sin(x) %输出精确解y1-y %输出偏差Adams法function Adams4 = Adams4()x = 0:pi/8:pi;h = pi/8;y = x;y(1) = 1;%% 由RK法求解Adams计算的初值for i = 1:3K1 = f(x(i),y(i));K2 = f(x(i)+0.5*h,y(i)+0.5*K1*h);K3 = f(x(i)+0.5*h,y(i)+0.5*K2*h);K4 = f(x(i)+h,y(i)+K3*h);y(i+1) = y(i)+h*(K1+2*K2+2*K3+K4)/6;endfor i = 4:length(x)-1s = 55*f(x(i),y(i))-59*f(x(i-1),y(i-1));s = s+37*f(x(i-2),y(i-2))-9*f(x(i-3),y(i-3)); %公式太长,分段写 y(i+1) = y(i)+s*h/24;s = 9*f(x(i+1),y(i+1))+19*f(x(i),y(i));s = s-5*f(x(i-1),y(i-1))+f(x(i-2),y(i-2));y(i+1) = y(i)+s*h/24;endy %输出Adams法的解y1 = cos(x)+sin(x) %输出精确解y1-y %输出偏差y的导函数function dy = f(a,b) dy = 2*cos(a)-b;。
数值分析习题与答案,DOC
第一章绪论习题一1.设x>0,x*的相对误差为δ,求f(x)=lnx的误差限。
解:求lnx的误差极限就是求f(x)=lnx的误差限,由公式即有2位有效数字,有5位有效数字,3.下列公式如何才比较准确?(1)(2)解:要使计算较准确,主要是避免两相近数相减,故应变换所给公式。
,利用:式计算误差最小。
四个选项:给定的数值表解:仍可使用n=1及n=2的Lagrange插值或Newton插值,并应用误差估计(5.8)。
线性插值时,用0.5及0.6两点,用Newton插值误差限,因,故二次插值时,用0.5,0.6,0.7三点,作二次Newton插值上给出的等距节点函数表,若用二次插值法求,函数表的步长令因得3.若,求和.解:由均差与导数关系于是4.若互异,求可知当有于是得求证.6.已知的函数表求出三次Newton均差插值多项式,计算f(0.23)的近似值并用均差的余项表达式估计误差.解:根据给定函数表构造均差表由于值并估计误差解:先构造差分表计算,用其中计算时用Newton后插公式(5.18)误差估计由公式(5.19)得这里仍为可先造,显然,再令令称为第二类求的表达式,并证明]上带权的正交解:因,即,故法方解得均方程为11.填空题(1)满足条件的插值多项式p(x)=().(2),则f[1,2,3,4]=(),f[1,2,3,4,5]=().(3)设为互异节点,为对应的四次插值基函数,则=(),=().是区间[其中,则==))))第4章数值积分与数值微分习题41.分别用复合梯形公式及复合Simpson公式计算下列积分.解本题只要根据复合梯形公式(6.11)及复合Simpson公式(6.13)直接计算即可。
对,取n=8,在分点处计算f(x)的值构造函数表。
按式(6.11)求出,按式(6.13)求得,积分)式估计误差,因,故(1)(2)(3)(1)令代入公式两端并使其相等,得解此方程组得,于是有再令,得故求积公式具有3次代数精确度。