清华大学高等数值分析 第一次实验作业

合集下载

清华大学数值分析A第一次作业

清华大学数值分析A第一次作业

7、设y0=28,按递推公式y n=y n−1−1100√783,n=1,2,…计算y100,若取√783≈27.982,试问计算y100将有多大误差?答:y100=y99−1100√783=y98−2100√783=⋯=y0−100100√783=28−√783若取√783≈27.982,则y100≈28−27.982=0.018,只有2位有效数字,y100的最大误差位0.00110、设f(x)=ln⁡(x−√x2−1),它等价于f(x)=−ln⁡(x+√x2−1)。

分别计算f(30),开方和对数取6位有效数字。

试问哪一个公式计算结果可靠?为什么?答:√x2−1≈29.9833则对于f(x)=ln(x−√x2−1),f(30)≈−4.09235对于f(x)=−ln(x+√x2−1),f(30)≈−4.09407而f(30)=⁡ln⁡(30−√302−1)⁡,约为−4.09407,则f(x)=−ln⁡(x+√x2−1)计算结果更可靠。

这是因为在公式f(x)=ln⁡(x−√x2−1)中,存在两相近数相减(x−√x2−1)的情况,导致算法数值不稳定。

11、求方程x2+62x+1=0的两个根,使它们具有四位有效数字。

答:x12=−62±√622−42=−31±√312−1则x1=−31−√312−1≈−31−30.98=−61.98x2=−31+√312−1=131+√312−1≈−131+30.98≈−0.0161312.(1)、计算√101.1−√101,要求具有4位有效数字答:√101.1−√101=√101.1+√101≈0.110.05+10.05≈0.00497514、试导出计算积分I n=∫x n4x+1dx1的一个递推公式,并讨论所得公式是否计算稳定。

答:I n=∫x n4x+1dx1 0=∫14(4x+1)x n−1−14x n−14x+1dx=114∫x n−11dx−14∫x n−14x+1dx1=14n−14I n−1,n=1,2…I0=∫14x+1dx=ln541记εn为I n的误差,则由递推公式可得εn=−14εn−1=⋯=(−14)nε0当n增大时,εn是减小的,故递推公式是计算稳定的。

数值分析大作业一

数值分析大作业一

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

JZX高等数值分析第一次实验作业

JZX高等数值分析第一次实验作业

相对残差 6.1302e-16 8.7797e-09 8.0295e-09 8.5677e-09 9.1433e-09
a、 m=1 (左为相对残差,右为取对数情况)
1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1
0 1
1.5
0
-5
-10
-15
-20
-25
-30
-35
-40
1000
4.761361
2.9675e-11
5.7069e-04
3.5336e+11
1500
15.646571
1.0778e-11
6.8236e-04
2.3672e+10
3000
132.198513
6.3164e-13
1.1700e-04
2.8110e+11
a、n=1000 时步数与相对残差关系图:(上为相对残差,下为取对数情况下结果)
(2)当 A 最大特征值远大于第二个特征值,最小特征值远小于第二个最小特征值时收敛
性情况。
思路:构造题目要求的矩阵 A。首先随机生成 n 阶矩阵 B,B 不满秩,构造对角阵 A1(最
大特征值远大于第二个最大特征值,最小特征值远小于第二个最小特征值),则由此构
造出对称正定矩阵 A: b1=B’*B; A=b1’*A1*b1。同样设定精确解 Xj 为元素全部为 1 的 n
5、 构造对称不定的矩阵,验证 Lanczos 方法的近似中断,观察收敛曲线中的峰点个数和特
征值的分布关系;观察当出现峰点时,MINRES 方法的收敛性态怎样。
解:思路:类似前两题,首先构造出一个 n 阶对角阵 D,其对角线上有 m 个负值,再对随

清华大学高等数值分析作业李津1——矩阵基础

清华大学高等数值分析作业李津1——矩阵基础

20130917题目求证:在矩阵的LU 分解中,111n n Tn ij i j j i j L I e e α-==+⎛⎫=- ⎪⎝⎭∑∑证明:在高斯消去过程中,假设0jj a ≠ ,若a=0,可以通过列变换使得前面的条件成立,这里不考虑这种情况。

对矩阵A 进行LU 分解,()()()()()1111111L M n M M M n ---=-=••-………… ,其中()1n Tn ij i j i j M j I e e α=+⎛⎫=+ ⎪⎝⎭∑ ,i e 、j e 为n 维线性空间的自然基。

()M j 是通过对单位阵进行初等变换得到,通过逆向的变换则可以得到单位阵,由此很容易得到()M j 的逆矩阵为1n T n ij i j i j I e e α=+⎛⎫- ⎪⎝⎭∑。

故111n n T n ij i j n j i j L I e e I α-==+⎛⎫⎛⎫=- ⎪ ⎪ ⎪⎝⎭⎝⎭∏∑上式中的每一项均是初等变换,从右向左乘,则每乘一次相当于对右边的矩阵进行一次向下乘法叠加的初等变换。

由于最初的矩阵为单位阵,变换从右向左展开,因而每一次变换不改变已经更新的数据,既该变换是从右向左一列一列更新数据,故11nn Tn ij i j j i j L I e e α==+⎛⎫=- ⎪⎝⎭∑∑。

数学证明:1n Tij i j i j e e α=+⎛⎫ ⎪⎝⎭∑具有,000n j j A -⎛⎫ ⎪⎝⎭ 和1,1000n j n j B -+-+⎛⎫⎪⎝⎭ 的形式,且有+1,-11,10000=000n j j n j n j A B --+-+⎛⎫⎛⎫⎪⎪⎝⎭⎝⎭而11n n T ij i j j k i j e e α-==+⎛⎫ ⎪⎝⎭∑∑具有1,1000n k n k B -+-+⎛⎫⎪⎝⎭的形式,因此: 1311111211121==n n n n n n T T T n ij i j n ij i j n ik i k j i j j i j k n i k n n T n i i n ik i i i k L I e e I e e I e e I e e I e ααααα---==+==+=-=+==+⎡⎤⎛⎫⎛⎫⎛⎫⎛⎫⎛⎫⎛⎫=---⎢⎥ ⎪ ⎪ ⎪ ⎪ ⎪⎪ ⎪ ⎪⎝⎭⎝⎭⎢⎥⎝⎭⎝⎭⎝⎭⎝⎭⎣⎦⎛⎫⎛⎫⎛⎫=-- ⎪ ⎪ ⎝⎭⎝⎝⎭∏∑∏∑∑∑∑∑……11211n n n T Tk n ik i kk k i k e I e e α--===+⎛⎫⎛⎫=- ⎪⎪ ⎪⎭⎝⎭⎝⎭∑∑∑#20130924题目一问:能否用逐次householder 相似变换变实矩阵A 为上三角矩阵,为什么?解:不能用逐次householder 相似变换变A 为上三角矩阵,原因如下:A 记作:()12=,,n A a a a ……, ,存在householder 阵1H s.t. 1111H a e α= ,则()()()111111111111111111111,,,0T Th H AH H a A H e H A H e H A H h H A H ααα⎛⎫'''=== ⎪⎪'⎝⎭⎛⎫''=+ ⎪ ⎪⎝⎭11H A H ''第一列的元素不能保证为1e 的倍数,故无法通过householder 变换实现上三角化。

数值分析实验报告-清华大学--线性代数方程组的数值解法

数值分析实验报告-清华大学--线性代数方程组的数值解法

数值分析实验报告-清华大学--线性代数方程组的数值解法(总15页)--本页仅作为文档封面,使用时请直接删除即可----内页可以根据需求调整合适字体及大小--线性代数方程组的数值解法实验1. 主元的选取与算法的稳定性问题提出:Gauss 消去法是我们在线性代数中已经熟悉的。

但由于计算机的数值运算是在一个有限的浮点数集合上进行的,如何才能确保Gauss 消去法作为数值算法的稳定性呢?Gauss 消去法从理论算法到数值算法,其关键是主元的选择。

主元的选择从数学理论上看起来平凡,它却是数值分析中十分典型的问题。

实验内容:考虑线性方程组 n n n R b R A b Ax ∈∈=⨯,,编制一个能自动选取主元,又能手动选取主元的求解线性方程组的Gauss 消去过程。

实验要求:(1)取矩阵⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡=1415157,6816816816 b A ,则方程有解T x )1,,1,1(* =。

取n=10计算矩阵的条件数。

让程序自动选取主元,结果如何?(2)现选择程序中手动选取主元的功能。

每步消去过程总选取按模最小或按模尽可能小的元素作为主元,观察并记录计算结果。

若每步消去过程总选取按模最大的元素作为主元,结果又如何?分析实验的结果。

(3)取矩阵阶数n=20或者更大,重复上述实验过程,观察记录并分析不同的问题及消去过程中选择不同的主元时计算结果的差异,说明主元素的选取在消去过程中的作用。

(4)选取其他你感兴趣的问题或者随机生成矩阵,计算其条件数。

重复上述实验,观察记录并分析实验结果。

程序清单n=input('矩阵A 的阶数:n=');A=6*diag(ones(1,n))+diag(ones(1,n-1),1)+8*diag(ones(1,n-1),-1); b=A*ones(n,1);p=input('计算条件数使用p-范数,p='); cond_A=cond(A,p) [m,n]=size(A);Ab=[A b];r=input('选主元方式(0:自动;1:手动),r=');Abfor i=1:n-1switch rcase(0)[aii,ip]=max(abs(Ab(i:n,i)));ip=ip+i-1;case (1)ip=input(['第',num2str(i),'步消元,请输入第',num2str(i),'列所选元素所处的行数:']);end;Ab([i ip],:)=Ab([ip i],:);aii=Ab(i,i);for k=i+1:nAb(k,i:n+1)=Ab(k,i:n+1)-(Ab(k,i)/aii)*Ab(i,i:n+1);end;if r==1Abendend;x=zeros(n,1);x(n)=Ab(n,n+1)/Ab(n,n);for i=n-1:-1:1x(i)=(Ab(i,n+1)-Ab(i,i+1:n)*x(i+1:n))/Ab(i,i);endx运行结果(1)n=10,矩阵的条件数及自动选主元Cond(A,1) =×103Cond(A,2) = ×103Cond(A,inf) =×103程序自动选择主元(列主元)a.输入数据矩阵A的阶数:n=10计算条件数使用p-范数,p=1选主元方式(0:自动;1:手动),r=0b.计算结果x=[1,1,1,1,1,1,1,1,1,1]T(2)n=10,手动选主元a. 每步消去过程总选取按模最小或按模尽可能小的元素作为主元矩阵A 的阶数:n=10计算条件数使用p-范数,p=1选主元方式(0:自动;1:手动),r=1(1)(1)61786115[]861158614A b ⎡⎤⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎣⎦第1步消元,请输入第1列所选元素所处的行数:1(2)(2) 6.0000 1.00007.00004.6667 1.0000 5.66678.0000 6.000015.0000[]8.00001.000015.00006.0000 1.00008.0000 6.0000 1.000015.00008.0000 6.000014.0000A b ⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦第2步消元,请输入第2列所选元素所处的行数:2…(实际选择时,第k 步选择主元处于第k 行) 最终计算得x=[, , , , , , , , , ]Tb. 每步消去过程总选取按模最大的元素作为主元 矩阵A 的阶数:n=10计算条件数使用p-范数,p=1选主元方式(0:自动;1:手动),r=1(1)(1)61786115[]861158614A b ⎡⎤⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎣⎦第1步消元,请输入第1列所选元素所处的行数:2(2)(2)8.0000 6.0000 1.000015.0000-3.50000.7500-4.250008.0000 6.0000 1.000015.0000[]8.0000 6.000015.00008.0000 1.00006.0000 1.000015.00008.0000 6.000014.0000A b ⎡⎤⎢⎥-⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦第2步消元,请输入第2列所选元素所处的行数:3…(实际选择时,第k 步选择主元处于第k+1行) 最终计算得x=[1,1,1,1,1,1,1,1,1,1]T(3)n=20,手动选主元a. 每步消去过程总选取按模最小或按模尽可能小的元素作为主元 矩阵A 的阶数:n=20计算条件数使用p-范数,p=1选主元方式(0:自动;1:手动),r=1(1)(1)61786115[]861158614A b ⎡⎤⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎣⎦第1步消元,请输入第1列所选元素所处的行数:1(2)(2) 6.0000 1.00007.00004.6667 1.0000 5.66678.0000 6.000015.0000[]8.00001.000015.00006.0000 1.00008.0000 6.0000 1.000015.00008.0000 6.000014.0000A b ⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦第2步消元,请输入第2列所选元素所处的行数:2…(实际选择时,第k 步选择主元处于第k 行) 最终计算得x=[,,,,,,,,,,,,,,,,,,,]T b. 每步消去过程总选取按模最大的元素作为主元 矩阵A 的阶数:n=20计算条件数使用p-范数,p=1选主元方式(0:自动;1:手动),r=1(1)(1)61786115[]861158614A b ⎡⎤⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎣⎦第1步消元,请输入第1列所选元素所处的行数:2(2)(2)8.0000 6.0000 1.000015.0000-3.50000.7500-4.250008.0000 6.0000 1.000015.0000[]8.0000 6.000015.00008.0000 1.00006.0000 1.000015.00008.0000 6.000014.0000A b ⎡⎤⎢⎥-⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦第2步消元,请输入第2列所选元素所处的行数:3…(实际选择时,第k步选择主元处于第k+1行)最终计算得x=[1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1]T(4)A分别为幻方矩阵,Hilbert矩阵,pascal矩阵和随机矩阵简要分析计算(1)表明:对于同一矩阵,不同范数定义的条件数是不同的;Gauss消去法在消去过程中选择模最大的主元能够得到比较精确的解。

《数值分析》课程设计—作业实验一...

《数值分析》课程设计—作业实验一...

《数值分析》课程设计—作业实验一1.1 水手、猴子和椰子问题:五个水手带了一只猴子来到南太平洋的一个荒岛上,发现那里有一大堆椰子。

由于旅途的颠簸,大家都很疲惫,很快就入睡了。

第一个水手醒来后,把椰子平分成五堆,将多余的一只给了猴子,他私藏了一堆后便又去睡了。

第二、第三、第四、第五个水手也陆续起来,和第一个水手一样,把椰子分成五堆,恰多一只猴子,私藏一堆,再去入睡,天亮以后,大家把余下的椰子重新等分成五堆,每人分一堆,正好余一只再给猴子,试问原先共有几只椰子?试分析椰子数目的变化规律,利用逆向递推的方法求解这一问题。

解:一、问题分析:对于本题,比较简单,我们只需要判断原来椰子的个数及每个人私藏了一份之后剩下的是否能被5除余1,直到最后分完。

二、问题求解:通过matlab 建立M 文件,有如下程序:或者对于第一个程序,n 取2000;对于第二个程序,n 取20001,就能得到我们想要的结果,即原先一共有15621个椰子,最终平均每人得4092个椰子。

n=input('input n:');forx=1:n p=5*x+1;for k=1:5 p=5*p/4+1;end if p==fix(p) break ; end enddisp([x,p]) input n:20001023 15621function fentao(n)a=cat(1,7);for j=n:-1:1 a(1)=j;i=1; while i<7a(i+1)=4*(a(i)-1)/5; i=i+1;endif a(7)==fix(a(7)) a, end end end>> fentao(20001) a =15621 12496 9996 7996 6396 5116 4092(本文档内的有些运行结果,限于篇幅,使文档结构更和谐、紧凑,已做相关的改动,程序代码没变)1.2 当0,1,2,,100n = 时,选择稳定的算法计算积分1d 10n xn xe I x e--=+⎰.解:一、问题分析:由1d 10n xn xe I x e--=+⎰知: 110101==+⎰dx I I以及: )1(110101011)1(1nnxxnxxn n n endx edx ee eI I ----+-+-==++=+⎰⎰得递推关系:⎪⎩⎪⎨⎧--=-=-+n n n I e n I I I 10)1(1101101,但是通过仔细观察就能知道上述递推公式每一步都将误差放大十倍,即使初始误差很小,但是误差的传播会逐步扩大,也就是说用它构造的算法是不稳定的,因此我们改进上述递推公式(算法)如下:⎪⎪⎩⎪⎪⎨⎧--=-=+-))1(1(101)1(101110n n n I e n I I I通过比较不难得出该误差是逐步缩小的,即算法是稳定的。

清华大学高等数值分析实验设计及答案

清华大学高等数值分析实验设计及答案

高等数值分析实验一工物研13 成彬彬2004310559一.用CG,Lanczos和MINRES方法求解大型稀疏对称正定矩阵Ax=b作实验中,A是利用A= sprandsym(S,[],rc,3)随机生成的一个对称正定阵,S是1043阶的一个稀疏阵A= sprandsym(S,[],0.01,3);检验所生成的矩阵A的特征如下:rank(A-A')=0 %即A=A’,A是对称的;rank(A)=1043 %A满秩cond(A)= 28.5908 %A是一个“好”阵1.CG方法利用CG方法解上面的线性方程组[x,flag,relres,iter,resvec] = pcg(A,b,1e-6,1043);结果如下:Iter=35,表示在35步时已经收敛到接近真实xrelres= norm(b-A*x)/norm(b)= 5.8907e-007为最终相对残差绘出A的特征值分布图和收敛曲线:S=svd(A); %绘制特征值分布subplot(211)plot(S);title('Distribution of A''s singular values');;xlabel('n')ylabel('singular values')subplot(212); %绘制收敛曲线semilogy(0:iter,resvec/norm(b),'-o');title('Convergence curve');xlabel('iteration number');ylabel('relative residual');得到如下图象:为了观察CG方法的收敛速度和A的特征值分布的关系,需要改变A的特征值:(1).研究A的最大最小特征值的变化对收敛速度的影响在A的构造过程中,通过改变A= sprandsym(S,[],rc,3)中的参数rc(1/rc为A的条件数),可以达到改变A的特征值分布的目的:通过改变rc=0.1,0.0001得到如下两幅图以上三种情况下,由收敛定理2.2.2计算得到的至多叠代次数分别为:48,14和486,由于上实验结果可以看出实际叠代次数都比上限值要小较多。

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

10
-10
0
100
200
300
400
500
600
700
800
900
迭代次数
图9
m=100时,Lanczos法求解Ax=b的收敛曲线
高等数值分析实验作业一
10
4
Lanzcos 算法的收敛曲线 (阶数 n=1002)
10
2
10
0
||rk||/||b||
10
-2
10
-4
10
-6
10
-8
10
-10
0
200
迭代次数
图12 m=10时,Minres法求解Ax=b的收敛曲线
10
2
Minres 算法的收敛曲线 (阶数 n=1002)
10
0
10
-2
||rk||/||b||
10
-4
10
-6
10
-8
10
-10
0
100
200
300
400
500
600
700
迭代次数
图13
10
2
m=50时,Minres法求解Ax=b的收敛曲线
10
0
Lanzcos 算法的收敛曲线 (阶数 n=1002)
m=10 m=50 m=100 m=400 m=800
10
-2
10
-4
||rk||/||b||
10
-6
10
-8
10
-10
10
-12
0
2
4
6
8
10
12
14
16
18
20
迭代次数
图6 对于b由m个特征向量线性表示时,Lanczos法求解Ax=b的收敛曲线
3. 当 A 只有 m 个不同特征值时, 对于大的 m 和小的 m, 观察有限精度下 Lanczos 方法如何收敛。 解:分别构建 m = 10、50、100、400、800 五个矩阵 A,分别求解 Ax=b,收敛 曲线如图 5 所示:
10
0
Lanzcos 算法的收敛曲线 (阶数 n=1002)
m=10 m=50 m=100 m=400 m=800
10
4
Lanzcos 算法的收敛曲线 (阶数 n=1002)
10
2
10
0
||rk||/||b||
10
-2
10
-4
10
-6
10
-8
0
50
100
150
200
250
300
350
400
迭代次数
图4
Lanczos法求解病态问题的收敛曲线
结论:对于良态正定问题, Lanczos法和CG法收敛性一样,相对残差的2范数随 迭代步数的增加而减小;对于病态正定问题 Lanczos法较CG法收敛的要慢一点, 而且收敛曲线更加不平滑,震荡较严重,相对残差的2范数同样没有最优性。 事实上,Lanczos法主要是解对称不定问题,针对这类问题才显现出它的优 越性,这也从侧面表明,对称正定问题不适合用Lanzcos法求解。
800
1000
1200
迭代次数
图16
m=500时,Minres法求解Ax=b的收敛曲线
结论:对于 Lanczos 方法,随着负特征值的增多,收敛曲线的峰点个数增多,振 荡越来越严重,发生近似中断的次数越来越多;然而,对于相同的 A,Minres 方法的相对残差没有出现峰值, 随着迭代数增加而单调下降。 收敛性态比 Lanczos 好,但代价是计算时间的急剧增大。
高等数值分析实验作业一
10
4
CG法收敛曲线 (阶数 n=1002)
10
2
10
0
||rk||/||b||
10
-2
10
-4
10
-6
10
-8
0
50
100
150
200
250
300
350
400
迭代次数
图2 CG法求解病态问题的收敛曲线
结论:对于良态问题,CG法收敛较快,相对残差的2范数随迭代步数的增加而减 小;对于病态问题CG法收敛的很慢,而且收敛曲线不平滑,震荡非常严重,相 对残差的2范数没有最优性。
2
10
0
||rk||/||b||
10
-2
10
-4
10
-6
10
-8
10
-10
0
100
200
300
400
500
600
700
迭代次数
图8
10
2
m=50时,Lanczos法求解Ax=b的收敛曲线
Lanzcos 算法的收敛曲线 (阶数 n=1002)
10
0
10
-2
||rk||/||b||
10
-4
10
-6
10
-8
10
0
CG法收敛曲线 (阶数 n=1002)
10
-1
10
-2
10
-3
||rk||/||b||
10
-4
10
-5
10
-6
10
-7
10
-8
0
20
40
60
80
100
120
140
160
180
迭代次数
图1 CG法求解良态问题的收敛曲线
(3)CG法求解Ax=b,A病态: 利用matlab编程实现CG算法。b = ones(1002,1),x0 = zeros(1002, 1)。计算 每一步迭代的残差rk相对于初始残差的2范数。相对残差2范数的对数值与 迭代步数的关系曲线如图2所示:
Minres 算法的收敛曲线 (阶数 n=1002)
10
0
10
-2
||rk||/||b||
10
-4
10
-6
10
-8
10
-10
0
100
200
300
400
500
600
700
800
900
迭代次数
图14
m=100时,Minres法求解Ax=b的收敛曲线
高等数值分析实验作业一
10
2
Minres 算法的收敛曲线 (阶数 n=1002)
结论:理论上,b 仅由 A 的 m 个不同特征向量的线性组合表示时,Lanczos 方法 必然 m 步收敛。但由于 A 的阶数为 1002,是比较大的,精度方面的限制导致计 算得到的 m 个特征向量并不都线性无关,所以,m = 800 时只需 20 步迭代。
高等数值分析实验作业一
5. 构造对称不定矩阵,验证 Lanczos 方法的近似中断,观察收敛曲线中的峰点 个数和特征值的分布关系; 观测当出现峰点时, MINRES 方法的收敛形态怎样。 解:分别构建负特征值个数为 m =10、50、100、200、500 的矩阵 A,分别计算 Ax = b 的解,Lanzcos 方法的收敛曲线如图 7-图 11 所示,Minres 的收敛曲线如 图 12-图 16 所示:
高等数值分析实验作业一
end
Lanczos 法 lanczos.m
function [T,Q,x,k,Errs,tiao]=Lanczos(A,b,x0,Err) x=[]; tiao=[]; [m,n]=size(b); if m<n b=b'; end [m,n]=size(x0); if m<n m=n; x0=x0'; end size(b) size(A*x0) r=b-A*x0; r_zeros=r; q=r/norm(r); q0=0; beta0=0; T=[0]; k=1; Q=[q]; Errs=[norm(r,2)/norm(b,2)]; while 1 disp('k=');disp(k); T=[T,zeros(k,1);zeros(1,k),0]; r=A*q-beta0*q0; T(k,k)=q'*r; r=r-T(k,k)*q; T(k+1,k)=norm(r,2); beta0=T(k+1,k); T(k,k+1)=beta0; q0=q; q=r/beta0; Tk=T(1:k,1:k) ; if k==1 disp('Ìø¹ý´Ë²½') tiao=[tiao,k]; else [L,U]=LanczosLU(Tk); bk=norm(r_zeros,2)*[1;zeros(k-1,1)];
高等数值分析实验作业一
附件:主要算法代码
CG 法 CG.m
function [x,Error,i,flag]=CG(A,b,x,ErrSet,uplimit) [m,n]=size(b); if m<n b=b'; end [m,n]=size(x); if m<n x=x'; end r=b-A*x; p=r; i=1; temp_rkrkplus=r'*r; Error=sqrt(temp_rkrkplus)/norm(b,2); while 1 temp_AP=A*p; temp_rkrk=temp_rkrkplus; temp_pAP=p'*temp_AP; if abs(temp_pAP)<1e-12 disp('¶ñÐÔÖжϣ¡') break; end a=temp_rkrk/(temp_pAP); x=x+a*p; r=r-a*temp_AP; temp_rkrkplus=r'*r; beta=temp_rkrkplus/temp_rkrk; p=r+beta*p; Err=sqrt(temp_rkrkplus)/norm(b,2); %/(norm(b)+temp_AP); if Err<ErrSet disp('Method converge£¡') disp(i) flag=1; break; end Error=[Error,sqrt(temp_rkrkplus)/norm(b,2)]; if i>uplimit flag=0; break end end end
相关文档
最新文档