数值分析作业

合集下载

数值分析大作业一

数值分析大作业一

数值分析大作业一一、算法设计方案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代表端点的一阶导。

数值分析第一章作业

数值分析第一章作业

习题(一)1. 指出四舍五入得到的下列各数有几位有效数字:x 1∗=7.8673,x 2∗=8.0916,x 3∗=0.06213,x 4∗=0.07800,x 5∗=90×103,x 6∗=2.0×10−4解:由有效数字定义得:x 1∗,x 2∗具有5位有效数字x 3∗,x 4∗具有4位有效数字x 5∗,x 6∗具有2位有效数字.2. 设准确值为x=3.78695,y=10,它们的近似值分别为x 1∗=3.7869,x 2∗=3.7870及y 1∗=9.9999,y 2∗=10.1,y 3∗=10.0001,试分析x 1∗,x 2∗,y 1∗,y 2∗,y 3∗分别具有几位有效数字. 解:x 1∗=3.7869=x 1∗=0.37869×101,k 1=1|x 1∗−x|=|3.7869−3.78695|=0.00005≤0.5×10−4=0.5×101−5, 即x 1∗具有5位有效数字;同理,x 2∗=3.7870=0.37870×101,k 2=1|x 2∗−x|=|3.7870−3.78695|=0.00005≤0.5×101−5,所以x 2∗具有5位有效数字; 将y 1∗,y 2∗,y 3∗分别写成y=±10k ×0.α1α2...αn 的表示形式,有:y 1∗=9.9999=0.99999×101,k 3=1;y 2∗=10.1=y 2∗=0.101×102,k 4=2;y 3∗=10.0001=0.100001×102,k 5=2;|y 1∗−y |=|9.9999−10|=0.0001=0.1×10−3≤0.5×101−4,n=4;|y 2∗−y |=|10.1−10|=0.1≤0.5×102−2,n=2;|y 3∗−y |=|10.0001−10|=0.0001=0.1×10−3≤0.5×102−5,n=5;所以y 1∗,y 2∗,y 3∗分别具有4,2,5位有效数字.8.为了使√11的近似值的相对误差不超过0.1%,问至少应取几位有效数字. 解:√11=0.3316624…=0.α1α2...αn ×10k ,α1=3,设x ∗有n 位有效数字,又因为|E x ∗|比值比较小, 故可用E r ∗(x ∗)= |E(x ∗)x ∗|代替相对误差E r ∗(x ∗),用εr ∗=εx ∗代替相对误差限εr 所以εr ∗≤12α1×10−n+1=16×10−n+1 令16×10−n+1≤0.1%,解得n ≥3.22即至少应取4位有效数字.12.如何计算下列函数值才比较精确.(1)11+2x −11+x ,对|x|≪1; (2)√x +1x −√x −1x ,对x ≫1;(3)∫dx 1+x 2N+1N,其中N 充分大; (4)1−cos xsin x ,对|x|≪1;(5)ln(30−√302−1)(开平方用6位函数表);解:(1)原式=1+x−(1+2x)(1+2x)(1+x)=−x (1+2x)(1+x); (2)原式=x+1x −(x−1x )√x+1x +√x−1x =2x √x+1x +√x−1x ;(3)原式=arc tan x|NN+1=arc tan N +1−arc tan N =arc tan N+1−N 1+N(N+1)=arc tan 11+N(N+1); (4)原式=2sinx 222sin x 2cos x 2=tan x2; (5)原式=30+√302−1=−ln(30+√302−1)令f(x)=ln(x −√x 2−1),则f(30)=ln(30−√302−1)=ln(30−√899),记a=30−√899 若用6位开方函数表,则有a ∗=30−29.9833=0.0167,故有ε(a ∗)=0.5×10−4, 而f(30)≈ln a ∗,于是ε(f (30))=ε(ln a ∗)≈|1a ∗|ε(a ∗)=0.50.0167×10−4≈0.003; 又因为f(x)等价于f(x)=-ln(x +√x 2−1),则f (30)=-ln(30+√899),记b=30+√899 同理b ∗=59.9833,进而ε(b ∗)=(2×10−4)−1,对f (30)≈ln b ∗ε(f (30))=ε(ln b ∗)≈|1b ∗|ε(b ∗)=0.559.9833×10−4≈0.834×10−6。

数值分析作业题(1)

数值分析作业题(1)

第一章 误差与算法1. 误差分为有__模型误差___, _观测误差___, __方法误差____, ___舍入误差____, Taylor 展开式近似表达函数产生的误差是_方法误差 .2. 插值余项是插值多项式的 方法误差。

0.2499作为1/4的近似值, 有几位有效数字?00.24990.249910,0m =⨯=即,031|0.2499|0.00010.5100.510,34m n n ---=<⨯=⨯=即22 3.1428751...,7=作为圆周率的近似值,误差和误差限分别是多少,有几位有效数字?2133.142875 3.14159260.00126450.5100.510---=<⨯=⨯有3位有效数字.* 有效数字与相对误差的关系3. 利用递推公式计算积分110,1,2,...,9n x n I x e dx n -==⎰错误!未找到引用源。

, 建立稳定的数值算法。

该算法是不稳定的。

因为:11()()...(1)!()n n n I n I n I εεε-=-==-111n n I I n n -=-, 10110I =4. 衡量算法优劣的指标有__时间复杂度,__空间复杂度_.时间复杂度是指: , 两个n 阶矩阵相乘的乘法次数是 , 则称两个n 阶矩阵相乘这一问题的时间复杂度为 .二 代数插值1.根据下表数据建立不超过二次的Lagrange 和Newton 插值多项式, 并写出误差估计式, 以及验证插值多项式的唯一性。

x 0 1 4f(x) 1 9 3Lagrange:设0120120,1,4;()1()9()3x x x f x f x f x ======则,, 对应 的标准基函数 为:1200102()()(1)(x 4)1()(1)(x 4)()()(01)(04)4x x x x x l x x x x x x ----===------ 1()...l x =2()...l x =因此, 所求插值多项式为:220()()()....i i i P x f x l x ===∑ (3)2()()(0)(1)(x 4)3!f R x x x ξ=--- Newton:构造出插商表:xi f(xi ) 一 二 三0 11 9 84 3 -2 -5/2所以, 所求插值多项式为:2001001201()()[,]()[,,]()()518(0)(0)(1)2...P x f x f x x x x f x x x x x x x x x x =+-+--=+----=插值余项: 2()[0,1,4,](0)(1)(x 4)R x f x x x =---2. 已知函数f(0)=1,f(1)=3,f(2)=7,则f[0,1]=___2________, f[0,1,2]=____1______)('],[000x f x x f =3.过0,1两节点构造三次Hermite 插值多项式, 使得满足插值条件: f(0)=1. .’(0)=... f(1.=2. .’(1)=1设0101010,1,()1()2'()0,'()1x x f x f x f x f x ======则,, 写出插商表:xi f(xi) 一 二 三0 10 1 01 a 1 11 a 1 0 a-1因此, 所求插值多项式为:插值余项:222()[0,0,1,1,](1)R x f x x x =-4.求f(x)=sinx 在[a,b]区间上的分段线性插值多项式, 并写出误差估计式。

数值分析大作业

数值分析大作业

数值分析大作业数值分析大作业姓名:黄晨晨学号:S1*******学院:储运与建筑工程学院学院班级:储建研17-2实验3.1 Gauss消去法的数值稳定性实验实验目的:理解高斯消元过程中出现小主元即很小时引起方程组解数值不定性实验内容:求解方程组Ax=b,其中(1)A1=0.3×10?1559.14315.291?6.130?1211.29521211,b1=59.1746.7812;(2)A2=10?7013 2.099999999999625?15?10102,b2=85.90000000000151;实验要求:(1)计算矩阵的条件数,判断系数矩阵是良态的还是病态的(2)用Gauss列主元消去法求得L和U及解向量x1,x2∈R4(3)用不选主元的高斯消去法求得L和U及解向量x1,x2∈R4(4)观察小主元并分析对计算结果的影响(1)计算矩阵的条件数,判断系数矩阵是良态的还是病态的代码:format longeformat compactA1=[0.3*10^-15,59.14,3,1;5.291,-6.130,-1,2;11.2,9,5,2;1,2,1,1] b1=[59.17;46.78;1;2]n=4C1=cond(A1,1) %C1为A1矩阵1范数下的条件数C2=cond(A1,2) %C2为A1矩阵2范数下的条件数C3=cond(A1,inf) %C3为1矩阵谱范数下的条件数结果:C1 =1.362944708720448e+02C2 =6.842955771253409e+01C3 =8.431146*********e+01显然A1矩阵为病态矩阵将矩阵A2,b2输入上述代码中求得A2矩阵的条件数为:C1 =1.928316831682894e+01C2 =8.993938090170119e+00C3 =1.835643564356072e+01显然A2矩阵也为病态矩阵(2)用Gauss列主元消去法求得L和U及解向量x1,x2∈R4代码:for k=1:n-1a=max(abs(A1(k:n,k)))if a==0returnendfor i=k:nif abs(A1(i,k))==ay=A1(i,:)A1(i,:)=A1(k,:)A1(k,:)=yx=b1(i,:)b1(i,:)=b1(k,:)b1(k,:)=xbreakendendif A1(k,k)~=0A1(k+1:n,k)=A1(k+1:n,k)/A1(k,k)A1(k+1:n,k+1:n)=A1(k+1:n,k+1:n)-A1(k+1:n,k)*A1(k,k+1:n) elsebreakendendL=tril(A1,0);for i=1:nL(i,i)=1;endLU=triu(A1,0)y1=L\b1x1=U\y1得到如下结果:x1 =3.845714853511634e+001.609517394778522e+00-1.547605454206655e+011.041130489899787e+01将A2=[10,-7,0,1;-3,2.0999********,6,2;5,-1,5,-1;0,1,0,2]b2=[8;5.900000000001;5;1]代入上述代码求得结果如下:X2 =4.440892098500626e-16-9.999999999999993e-019.999999999999997e-011.000000000000000e+00(3)用不选主元的高斯消去法求得L和U及解向量x1,x2∈R4代码:format longeformat compactA1=[0.3*10^-15,59.14,3,1;5.291,-6.130,-1,2;11.2,9,5,2;1,2,1,1] b1=[59.17;46.78;1;2][L,U]=lu(A1)y1=L\b1x1=U\y1求得如下结果:x1=3.845714853511634e+001.609517394778522e+00-1.547605454206655e+011.041130489899787e+01将A2=[10,-7,0,1;-3,2.0999********,6,2;5,-1,5,-1;0,1,0,2] b2=[8;5.900000000001;5;1]代入上述代码,求得结果如下:x 2 =4.440892098500626e-16 -9.999999999999993e-01 9.999999999999997e-01 9.999999999999999e-01(2)(3)求得结果相同,可验证结果正确。

数值分析作业

数值分析作业

数值分析作业及答案Chap11、写出下列语句的运行结果。

在MA TLAB 上执行它们以验证所得解答。

a=[1 2 3 ;4 5 6 ]’ b=[9;7;5;3;1] c=b(2:4) d=b(4:-1:1) e=sort(b) f=[3,b ’]解:a=635241 b=13579c=357d=9753 e=97531F=[3 9 7 5 3 1] 3、给定一向量:a=[4 -1 2 -8 4 5 -3 -1 6 -7]写一段程序计算a 中正数的和。

运行程序并显示结果。

解:a=[4 -1 2 -8 4 5 -3 -1 6 -7]; s=0;for i=1:length(a) if a(i)>0s=s+a(i); end end s6、编写一个函数M 文件fun_es(x),计算如下函数:230.5sin x y e x x =-其中参数可以为标量,也可以为向量。

在MA TLAB 里键入如下命令检验此函数:fun_es(3) fun_es([1 2 3])解:function y=fun_es(x) y=0.5*exp(x/3)-x.^2.*sin(x);chap21、设0x >,x 的相对误差为δ,求L nx 的误差。

解:Lnx-Lnx*=dLnx=dx/x=δ2、设x 的相对误差为2%,求2x 的相对误差。

解:dLnf(x)=xf ’(x)/f(x)dLnx=4%5、计算球体积要使相对误差限为1%,问度量半径R 时允许的相对误差限是多少?解:dLnf(x)=xf ’(x)/f(x)dLnx=3dLnx=1% dLnx=0.33%9、正方形的边长大约为100cm ,应怎样测量才能使其面积误差不超过1cm 2? 解:s=x 2s-s*=2x(x-x*)=1x-x*=1/(2x)=1/200=0.5*10-2 即测量边的误差不超过0.005cm 10、设212S gt =,假定g 是准确的,而对t 的测量有±0.1秒的误差,证明当t 增加时S 的绝对误差增加,而相对误差却减少。

(完整版)数值分析第一次作业

(完整版)数值分析第一次作业

问题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代表端点的一阶导。

数值分析练习题附答案

数值分析练习题附答案

目录一、绪论------------------------------------------------------------------------------------- 2-2二、线性方程组直接解法列主元高斯LU LDL T GG T-------------------- 3-6二、线性方程组迭代法----------------------------------------------------------------- 7-10 三、四、非线性方程组数值解法二分法不动点迭代---------------------- 11-13五、非线性方程组数值解法牛顿迭代下山弦截法----------------- 14-15六、插值线性插值抛物线插值------------------------------------------------ 16-18七、插值Hermite插值分段线性插值-----------------------------------------19-22八、拟合------------------------------------------------------------------------------------ 23-24九、数值积分----------------------------------------------------------------------------- 25-29十、常微分方程数值解法梯形欧拉改进----------------------------------- 30-32 十一、常微分方程数值解法龙格库塔------------------------------------------ 33-35绪论1-1 下列各数都是经过四舍五入得到的近似值 ,试分别指出它们的绝对误差限,相对误差限和有效数字的位数.X 1 =5.420, X 2 =0.5420, X 3 =0.00542, X 4 =6000, X 5 =0.6×105注:将近似值改写为标准形式X 1 =(5*10-1+4*10-2+2*10-3+0*10-4)*101 即n=4,m=1 绝对误差限|△X 1|=|X *1-X 1|≤ 12×10m-n =12×10-3 相对误差限|△r X 1|= |X∗1−X1||X∗1|≤|X∗1−X1||X1|= 12×10-3/5.4201-2 为了使101/2 的相对误差小于0.01%, 试问应取几位有效数字?1-3 求方程x 2 -56x+1=0的两个根, 使它们至少具有4位有效数字( √783≈27.982)注:原方程可改写为(x-28)2=783线性方程组解法(直接法)2-1用列主元Gauss消元法解方程组解:回代得解:X1=0 X2=-1 X3=12-2对矩阵A进行LU分解,并求解方程组Ax=b,其中解:(注:详细分解请看课本P25)A=(211132122)→(211(1/2)5/23/2(1/2)3/23/2)→(2111/25/23/21/2(3/5)3/5)即A=L×U=(11/211/23/51)×(2115/23/23/5)先用前代法解L y=P b 其中P为单位阵(原因是A矩阵未进行行变换)即L y=P b 等价为(11/211/23/51)(y1y2y3)=(111)(465)解得 y 1=4 y 2=4 y 3=35再用回代解Ux =y ,得到结果x即Ux =y 等价为(2115/23/23/5)(x 1x 2x 3)=(y 1y 2y 3)=(443/5) 解得 x 1=1 x 2=1 x 3=1即方程组Ax=b 的解为x =(111)2-3 对矩阵A 进行LDL T 分解和GG T 分解,求解方程组Ax=b,其中A=(164845−48−422) , b =(123)解:(注:课本 P 26 P 27 根平方法)设L=(l i j ),D=diag(d i ),对k=1,2,…,n,其中d k =a kk -∑l kj 2k−1j=1d jl ik =(a ik −∑l ij l kj k−1j=1d j )/ d k 即d 1=a 11-∑l 1j 20j=1d j =16-0=16因为 l 21=(a 21−∑l 2j l 1j 0j=1d j )/ d 1=a 21/ d 1=416=14 所以d 2=a 22-∑l 2j 21j=1d j =5-(14)2d 1=4同理可得d 3=9 即得 D=(1649)同理l 11=(a 11−∑l ij l 1j 0j=1d j )/ d 1=1616=1=l 22=l 33 l 21=(a 21−∑l 2j l 1j 0j=1d j )/ d 1=416=14 l 31=(a 31−∑l 3j l 1j 0j=1d j )/ d 1=816=12 l 32=(a 32−∑l 3j l 2j 1j=1d j )/ d 2=−4−12×14×164=−64=-32即L=(114112−321) L T=(114121−321) 即LDL T分解为A=(114112−321)(1649)(114121−321)解解:A=(164845−48−422)→(41212−32−33)故得GG T分解:A=(4122−33)(4122−33) LDL T分解为A=(114112−321)(1649)(114121−321) 由(114112−321)(y 1y 2y 3)=(123) ,得(y 1y 2y 3)=(0.250.8751.7083)再由(4122−33)(x 1x 2x 3)=(0.250.8751.7083) ,得(x 1x 2x 3)=(−0.54511.29160.5694)2-4 用追赶法求解方程组:解:(4−1−14−1−14−1−14−1−14)→(4−14−1154−415−15615−1556−120956−56209−1780209)由(4−1154−15615−120956−1780209)(y1y2y3y4y5)=(100200),得(y1y2y3y4y5)=(256.66671.785700.4784753.718)再由(1−141−4151−15561−562091)(x1x2x3x4x5)=(256.66671.785700.4784753.718),得(x1x2x3x4x5)=(27.0518.20525.769314.87253.718)线性方程组解法(迭代法)2-1 设线性方程组{4x 1−x 2+2x 3=1−x 1−5x 2+x 3=22x 1+x 2+6x 3=3(1) 写出Jacobi 法和SOR 法的迭代格式(分量形式) (2) 讨论这两种迭代法的收敛性(3) 取初值x (0)=(0,0,0)T ,若用Jacobi 迭代法计算时,预估误差 ||x*-x (10)||∞ (取三位有效数字)解:(1)Jacobi 法和SOR 法的迭代格式分别为Jacobi 法迭代格式SOR(2)因为A 是严格对角占优矩阵,但不是正定矩阵,故Jacobi 法收敛,SOR 法当0<ω≤1时收敛.⎪⎪⎪⎩⎪⎪⎪⎨⎧+--=-+-=+-=+++216131525151412141)(2)(1)1(3)(3)(1)1(2)(3)(2)1(1k k k k k k k k k x x x x x x xx x ⎪⎪⎪⎩⎪⎪⎪⎨⎧-++-=+-+-=+-+-+=++++++)216131()525151()412141()(3)1(2)1(1)(3)1(3)(3)(2)1(1)(2)1(2)(3)(2)(1)(1)1(1k k k k k k k k k k k k k k k x x x x x x x x x x x x x x x ωωω(3)由(1)可见||B ||∞=3/4,且取x (0)=(0,0,0)T ,经计算可得x (1)=(1/4,-2/5,1/2)T ,于是||x (1)-x (0)||∞=1/2,所以有2-2 设方程组为{5x 1+2x 2+x 3=−12−x 1+4x 2+2x 3=202x 1−3x 2+10x 3=3试写出其Jacobi 分量迭代格式以及相应的迭代矩阵,并求解。

数值分析(在线作业)

数值分析(在线作业)

数值分析(在线作业)单选:1、设是经过四舍五入后得到的近似值,则分别有几位有效数字?(A )A、3,3B、2,4C、3,4D、4,32、计算球的体积时,为使其相对误差限为1%,测量半径R时,相对误差最大为(B )A、1%B、0.33%C、3%D、9.9%3、超定方程组的最小二乘解为(A)A、B、C、D、4、已知则为(D)A、2B、6C、-6D、85、已知A=则为(D)A、0.367B、0C、-34D、396016、设方程组Ax=b,其中则A能进行Cholesky分解(即A=LLT,其中L为下三角矩阵)时,取值范围为(A )A、B、C、D、7、设,则差商为(A )A、1B、-1C、0D、28、设,则差商为(C)A、1B、-1C、0D、29、已知函数表为分别用Newton向前、向后插值公式计算f(1.5),f(3.7)的近似值(B )A、B、C、D、10、设,则的Newton迭代公式为(A )A、B、C、D、11、设,则当的Newton迭代收敛时,的取值范围为(A)A、B、C、D、12、已知325413有6位有效数字,则绝对误差限为(B )A、0.05B、0.5C、0.005D、513、已知,则下列哪个多项式为的二次最佳平方逼近(B)A、B、C、D、14、计算积分,若用复合Simpson公式进行近似计算,并且想误差不超过则至少要进行多少等分?( C )A、3B、6C、12D、2415、给定线性方程组,其中,,使用迭代公式,若迭代收敛,则的取值范围为(D)A、B、C、D、16、已知,则的谱半径为(C )A、3B、1C、7D、817、已知,则的拉格朗日插值多项式为(A )A、B、C、D、18、已知,则用梯形公式计算积分与精确值相比(A)A、偏大B、偏小C、相等D、不确定19、对于线性方程组,则雅克比迭代与高斯-赛德尔迭代的敛散性分别为(A )A、收敛发散B、收敛收敛C、发散收敛D、发散发散20、已知函数值,则均差为(D )A、B、6 C、10 D、2。

数值分析最佳习题(含答案)

数值分析最佳习题(含答案)

第一章 绪论姓名 学号 班级习题主要考察点:有效数字的计算、计算方法的比较选择、误差和误差限的计算。

1 若误差限为5105.0-⨯,那么近似数0.003400有几位有效数字?(有效数字的计算) 解:2*103400.0-⨯=x ,325*10211021---⨯=⨯≤-x x 故具有3位有效数字。

2 14159.3=π具有4位有效数字的近似值是多少?(有效数字的计算) 解:10314159.0⨯= π,欲使其近似值*π具有4位有效数字,必需41*1021-⨯≤-ππ,3*310211021--⨯+≤≤⨯-πππ,即14209.314109.3*≤≤π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 =的相对误差。

北航研究生数值分析编程大作业1

北航研究生数值分析编程大作业1

数值分析大作业一、算法设计方案1、矩阵初始化矩阵[]501501⨯=ij a A 的下半带宽r=2,上半带宽s=2,设置矩阵[][]5011++s r C ,在矩阵C 中检索矩阵A 中的带内元素ij a 的方法是:j s j i ij c a ,1++-=。

这样所需要的存储单元数大大减少,从而极大提高了运算效率。

2、利用幂法求出5011λλ,幂法迭代格式:0111111nk k k k kk T k k k u R y u u Ay y u ηηβ------⎧∈⎪⎪=⎪=⎨⎪=⎪⎪=⎩非零向量 当1210/-≤-k k βββ时,迭代终止。

首先对于矩阵A 利用幂法迭代求出一个λ,然后求出矩阵B ,其中I A B λ-=(I 为单位矩阵),对矩阵B 进行幂法迭代,求出λ',之后令λλλ+'='',比较的大小与λλ'',大者为501λ,小者为1λ。

3、利用反幂法求出ik s λλ,反幂法迭代格式:0111111nk k k k kk T k k k u R y u Au y y u ηηβ------⎧∈⎪⎪=⎪=⎨⎪=⎪⎪=⎩非零向量 当1210/-≤-k k βββ时,迭代终止,1s k λβ=。

每迭代一次都要求解一次线性方程组1-=k k y Au ,求解过程为:(1)作分解LU A =对于n k ,...,2,1=执行[][]s k n r k k k i c c c c c n s k k k j c cc c k s ks k t k s k r i t t s t i k s k i k s k i js j t k s j r k t t s t k j s j k j s j k <+++=-=++=-=+++----=++-++-++-++----=++-++-++-∑∑);,min(,...,2,1/)(:),min(,...,1,:,1,11),,1max(,1,1,1,11),,1max(,1,1,1(2)求解y Ux b Ly ==,(数组b 先是存放原方程组右端向量,后来存放中间向量y))1,...,2,1(/)(:/:),...,3,2(:,1),min(1.1.11),1max(,1--=-===-=+++-++-+--=++-∑∑n n i c x c b x c b x n i b c b b i s t n s i i t t s t i i i ns n n ti r i t t s t i i i使用反幂法,直接可以求得矩阵按模最小的特征值s λ。

北理工数值分析大作业

北理工数值分析大作业

数值分析上机作业第 1 章1.1计算积分,n=9。

(要求计算结果具有6位有效数字)程序:n=1:19;I=zeros(1,19);I(19)=1/2*((exp(-1)/20)+(1/20));I(18)=1/2*((exp(-1)/19)+(1/19));for i=2:10I(19-i)=1/(20-i)*(1-I(20-i));endformat longdisp(I(1:19))结果截图及分析:在MATLAB中运行以上代码,得到结果如下图所示:当计算到数列的第10项时,所得的结果即为n=9时的准确积分值。

取6位有效数字可得.1.2分别将区间[-10.10]分为100,200,400等份,利用mesh或surf命令画出二元函数z=的三维图形。

程序:>> x = -10:0.1:10;y = -10:0.1:10;[X,Y] = meshgrid(x,y);Z = exp(-abs(X))+cos(X+Y)+1./(X.^2+Y.^2+1);subplot(2,2,1);mesh(X,Y,Z);title('步长0.1')>> x = -10:0.2:10;y = -10:0.2:10;[X,Y] = meshgrid(x,y);Z = exp(-abs(X))+cos(X+Y)+1./(X.^2+Y.^2+1);subplot(2,2,1);mesh(X,Y,Z);title('步长 0.2')>>x = -10:0.05:10;y = -10:0.05:10;[X,Y] = meshgrid(x,y);Z = exp(-abs(X))+cos(X+Y)+1./(X.^2+Y.^2+1);subplot(2,2,1);mesh(X,Y,Z);title('步长0.05')结果截图及分析:由图可知,步长越小时,绘得的图形越精确。

数值分析期末大作业

数值分析期末大作业

一、问题提出设方程f(x)=x 3-3x-1=0有三个实根 x *1=1.8793 , x *2=-0.34727 ,x *3=-1.53209现采用下面六种不同计算格式,求 f(x)=0的根 x *1 或x *2 。

1、 x = 213xx + 2、x = 313-x3、 x = 313+x4、 x = 312-x 5、 x = x13+6、 x = x - ()1133123---x x x二、目的和意义1、通过实验进一步了解方程求根的算法;2、认识选择计算格式的重要性;3、掌握迭代算法和精度控制;4、明确迭代收敛性与初值选取的关系。

三、结构程序设计本程序实在matlab 软件上进行操作的。

首先建立一个空白的M-文件。

在编辑器中输入以下内容,并保存。

function [X1,m,n,q]=shizi1(p) x=zeros(100,1); x=double(x);x(1,1)=p;i=1;deltax=100;while (i<100 & deltax > 0.000001)x(i+1,1)=(3*x(i,1)+1)/x(i,1)^2deltax=abs(x(i+1,1)-x(i,1));i=i+1;endX1=x(1,1);m=i;n=x(i,1);q=deltax;以上是运行函数,下一步在建立一个执行M-文件,输入以下内容,并保存。

其中X1为初始值,m为迭代次数,n为最后得到的值,q为|x k+1-x k|。

clear all;clc;p=1.8;[X1,m,n,q]=shizi1(p)1、对第一个迭代公式,在执行文件中输入p=1.8;[X1,m,n,q]=shizi1(p)。

得到如下结果如下:初值为1.8,迭代100次,精度为10-6。

可见该迭代公式是发散的,将初值改为-1.5,其他均条件不变。

p=-1.5;[X1,m,n,q]=shizi1(p)改变初值后可以得到一个接近真值的结果x*3的结果ans=-1.5321。

数值分析习题(含标准答案)

数值分析习题(含标准答案)

数值分析习题(含标准答案)
一、选择题(每题5分,共20分)
1. 下列哪个选项不属于数值分析的研究范畴?
A. 数值微分
B. 数值积分
C. 数值逼近
D. 数据库管理
答案:D
2. 在数值分析中,求解线性方程组常用的方法有?
A. 高斯消元法
B. 迭代法
C. 拉格朗日乘数法
D. 上述所有方法
答案:D
3. 下列哪种方法适用于求解非线性方程组?
A. 牛顿法
B. 梯度下降法
C. 高斯消元法
D. 上述所有方法
答案:D
4. 在数值积分中,下列哪种方法具有最高的精度?
A. 梯形法则
B. 辛普森法则
C. 高斯求积法
D. 上述所有方法
答案:C
二、填空题(每题5分,共20分)
1. 数值分析的主要目的是通过有限步骤的运算,对数学问题进行近似求解。

2. 在数值微分中,常用的差分公式有前向差分、后向差分和中心差分。

3. 数值逼近的主要方法包括插值法和逼近法。

4. 在数值积分中,常用的方法有梯形法则、辛普森法则和高斯求积法。

三、解答题(每题10分,共30分)
1. 已知函数 f(x) = e^x,求其在 x = 0.5 处的导数。

答案:f'(0.5) ≈ 1.6487
2. 求解线性方程组 2x + 3y = 5,4x y = 1。

答案:x ≈ 0.625,y ≈ 1.25
3. 已知函数 f(x) = x^3 3x^2 + 4,求其在区间 [0, 2] 上的积分。

答案:f(x) 在区间 [0, 2] 上的积分≈ 3.6667。

数值分析上机作业(2)

数值分析上机作业(2)

一、数值求解如下正方形域上的Poisson 方程边值问题 2222(,)1,0,1(0,)(1,)(1),01(,0)(,1)0,01u u f x y x y x y u y u y y y y u x u x x ⎧⎛⎫∂∂-+==<<⎪ ⎪∂∂⎪⎝⎭⎨==-≤≤⎪⎪==≤≤⎩二、用椭圆型第一边值问题的五点差分格式得到线性方程组为2,1,1,,1,10,1,,0,141,?,?,?,?0,1i j i j i j i j i j ijj N j i i N u u u u u h f i j N u u u u i j N -+-+++----=≤≤====≤≤+, 写成矩阵形式Au=f 。

其中1.三 、编写求解线性方程组Au=f 的算法程序, 用下列方法编程计算, 并比较计算速度。

2.用Jacobi 迭代法求解线性方程组Au=f 。

3.用块Jacobi 迭代法求解线性方程组Au=f 。

4. 用SOR 迭代法求解线性方程组Au=f,用试算法确定最佳松弛因子。

1122N N v b v b u f v b ⎛⎫⎛⎫ ⎪ ⎪ ⎪ ⎪== ⎪ ⎪ ⎪ ⎪⎝⎭⎝⎭4114114ii A -⎛⎫ ⎪- ⎪= ⎪- ⎪-⎝⎭11,12,1,121,22,2,21,2,,2211,12,1,121,22,2,221,2,,(,,...,),(,,...,),......,(,,...,)(,,...,)?,(,,...,)?,......,(,,...,)?1,999,0.10.011T T N N TN N N N N T T N N T N N N N N v u u u v u u u v u u u b h f f f b h f f f b h f f f h N h N ====+=+=+===+取或则或,1,,1,2,...,i j f i j N== 1122NN A I I A A I I A -⎛⎫ ⎪- ⎪= ⎪- ⎪-⎝⎭5.用块SOR 迭代法求解线性方程组Au=f,用试算法确定最佳松弛因子。

数值分析作业答案——第五版

数值分析作业答案——第五版

数值分析第一次作业及参考答案1. 已测得函数()y f x =的三对数据:(0,1),(-1,5),(2,-1),(1)用Lagrange 插值求二次插值多项式。

(2)构造差商表。

(3)用Newton 插值求二次插值多项式。

解:(1)Lagrange 插值基函数为0(1)(2)1()(1)(2)(01)(02)2x x l x x x +-==-+-+-同理 1211()(2),()(1)36l x x x l x x x =-=+故2202151()()(1)(2)(2)(1)23631i i i p x y l x x x x x x x x x =-==-+-+-++=-+∑(2)令0120,1,2x x x ==-=,则一阶差商、二阶差商为0112155(1)[,]4,[,]20(1)12f x x f x x ---==-==-----0124(2)[,,]102f x x x ---==-22()1(4)(0)1*(0)(1)31P x x x x x x =+--+-+=-+2. 在44x -≤≤上给出()x f x e =的等距节点函数表,若用二次插值求xe 的近似值,要使截断误差不超过610-,问使用函数表的步长h 应取多少?解:()40000(),(),[4,4],,,, 1.xk x f x e fx e e x x h x x h x x th t ==≤∈--+=+≤考察点及(3)200044343()()[(()]()[()]3!(1)(1)(1)(1)3!3!2.(4,4).6fR x x x h x x x x ht t tet h th t h e heξξ=----+-+≤+⋅⋅-=≤∈-则436((1)(1)100.006.t t th h--+±<<在点取到极大值令 得3.求2()f x x=在[a,b]上的分段线性插值函数()hI x,并估计误差。

  1. 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
  2. 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
  3. 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。

数值分析课后作业:习题一
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.下列各近似值的绝对误差限都是102
1⨯-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+b
ans =
5 7 9
>> b-a
ans =
3 3 3
>> a.*b
ans =
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];
>> c
c = 1 2 3 1 2 3 1 2 3
d = 4 5 6 4 5 6 4 5 6 >> c+d
ans =
5 7 9 5 7 9 5 7 9
>> d-c
ans = 3 3 3
3 3 3
3 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=1
ax c
1 x>0 其中 sgn = 0 x=0 -1 x<0 disp('Please input the coefficients of');
disp('quadratic equation ax^2+bx+c=0, resp
ectively') 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 End
delta=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>0
x 1=(-b-temp)/(2*a) End if b<0
x 1=(-b+temp)/(2*a) End if b=0
x 1=temp/(2*a) End
x 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)<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 End
delta=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>0
x 1=digit(digit(-b-temp,m)/digit(2*a,m),m) ; End if b<0
x 1=digit(digit(-b+temp,m)/digit(2*a,m),m); End if b=0
x 1=digit(temp/digit(2*a,m),m); End
x 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(1
1
≤<--+∞
=∑x n
x n
n n 和ln
11...),1
2...53(2111
253<<-++++++=-++x n x x x x x x n ,给出计算ln2的近似方法,编写相应的Matlab 程序,并比较算法运行情况. 解:方法一: x=1; s=0;
for k=1:100
s=s+(-1)^(k+1)*(x^k)/k; end s
q=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
所以方法二较方法一好。

相关文档
最新文档