数值分析计算实习题一
学号：
姓名：
院系：
2015年11月5日
一、分析
1.1算法分析
题目要求求出：
1）特征值从小到大排列的最小特征值1λ和最大特征值501λ。 2）特征值中模最小的特征值s λ。 3）靠近一组数k μ的一组特征值k i λ。 4）矩阵A 的条件数cond(A)2。 5）行列式detA 。 解决方法：
1）若将所有行列式按模的大小排列则模最大的特征值一定是1λ和501λ中的一个，因此利用幂法求出模最大的特征值1m λ。然后利用带原点平移的幂法，将系数矩阵变为1m A I λ-即将所有特征值都减去1m λ，则特征值按大小顺序排列的次序不变，模最大的特征值依然在整个排列的两端，再用一次幂法得到模最大的特征值21=m m λλλ-，
其中λ为带原点平移的幂法求出的特征值，最后两个特征值1m λ、2m λ比较大小，大的为501λ，小的为1λ。
2）因为s λ为按模最小的特征值，因此用反幂法可求的其特征值。
3）因为k i λ靠近数k μ，因此k i k λμ-一定是所有的k λμ-中模最小的，因此可利用带原点平移的反幂法求出特征值k i λ，此时的系数矩阵变为k A I μ-。
4）条件数cond(A)2为模最小的特征值与模最大的特征值的比的绝对值，因此利用1和2中求出的1m λ和s λ可解出条件数。
5）可对矩阵A 进行LU 分解，即A LU =则det()det()det()A L U =?，又因为矩阵
L 对角线元素为1，则det()L
=1，所以det()det()A U =，U 为上三角阵，行列式为对角线元素的乘积，因此可得A 的行列式。
1.2程序分析
1.2.1 因为A 为拟三角阵，储存时零元素不储存，因此将矩阵A 压缩为5*501的矩阵C
A 的带内元素ij a =C 中的元素1,i j s j c -++ 程序中A[5][501]即为压缩后的矩阵。
1.1.2 程序中的B[5][501]为过渡矩阵，在幂法迭代、反幂法迭代以及LU 分解中均用矩阵B 来计算，计算之间对B 进行适当的赋值。先令B=A ，调用幂法函数求1m λ；再令B=A-1m λI ，调用幂法函数求出2m λ；再令B=A ，调用反幂法函数求
s λ；在令B=k A I μ-，调用反幂法函数求k
i λ，最后令B=A ，将B 进行LU 分解，
计算A 的行列式。 1.2.3幂法求解过程： 1）取初始向量u 0=[1,1,…1]; 2）进行k 次迭代，k=1,2,3…
3）对每一次迭代，计算上一次迭代中的(1)k u -的模（2范数）1k η-=将(1)k u -单位化后赋值给1k y -，即1(1
)1/k k k
y u η
---
=，计算矩阵B 与向量1k y -的乘积，
1k k u By -=，计算1k y -与k u 的内积1T k k k y u β-=。
4）如果两次相邻两次迭代的β满足：1/k k k βββε--≤（允许误差），则结束迭代，k β的值可认为是矩阵B 对应的模最大的特征值。如果不满足误差条件，则重复3）的迭代直到达到误差允许值。 1.2.4 反幂法求解过程：
1）取初始向量u 0=[1,1,…1];对矩阵B 进行LU 分解。 2）进行k 次迭代，k=1,2,3…
3）对每一次迭代，计算上一次迭代中的(1)k u -的模（2范数）1k η-=将(1)k u -单位化后赋值给1k y -，即1(1)1/k k k y u η---=，利用LU 分解后的矩阵B ，求
线性方程组1k k Bu y -=，计算1k y -与k u 的内积1T k k k y u β-=。
4）如果两次相邻两次迭代的β满足：1/k k k βββε--≤（允许误差），则结束迭
的值可认为是矩阵B对应的模最大的特征值。如果不满足误差条件，则代，
k
重复3）的迭代直到达到误差允许值。
1.2.5 拟三角矩阵的LU分解过程可参照数值分析书26页的算法得到。
二、程序
整个数值分析在c语言中的程序如下：
#include
#include
doublemifa();
doublefanmifa();
void fenjie();//声明三个函数，mifa是用幂法求模最大的特征值，fanmifa是用反幂法求模最小的特征值，fenjie是对系数矩阵做LU分解
double B[5][501]={0};
double A[5][501]={0};//定义两个系数矩阵为全局变量,B作为求解过程中的系数矩阵使用
double sigma=1e-12;//定义误差允许值
/********主程序********/
void main()
{
inti,j,t;
doublelamk,lamm,min,max,miu,lam,x,lams;
double cond2;
doubledet=1;
for(i=0;i<=500;i++)//输入矩阵A
{
j=i+1;
x=0.1/j;
A[2][i]=(1.64-0.024*j)*sin(0.2*j)-0.64*exp(x);
A[3][i]=0.16;
A[4][i]=-0.064;
}
for(i=1;i<=500;i++)
{
A[1][i]=0.16;
}
for(i=2;i<=500;i++)
{
A[0][i]=-0.064;
}
A[3][500]=0;
A[4][500]=0;
A[4][499]=0;
for(i=0;i<=4;i++)//先使B=A
{
for(j=0;j<=500;j++)
{
B[i][j]=A[i][j];
}
}
lamk=mifa();//幂法求模最大的特征值
for(i=0;i<=4;i++)//赋值计算B=（A-lamk*I）
{
for(j=0;j<=500;j++)
{
if(i==2) B[i][j]=A[i][j]-lamk;
else B[i][j]=A[i][j];
}
}
lamm=mifa()+lamk;//带原点位移的幂法求位移后模最大的特征值if(lamk else {min=lamm;max=lamk;} printf("lam1=%13.11le\n",min);//显示12位有效数字 printf("lam501=%13.11le\n",max); for(i=0;i<=4;i++)//赋值B=A { for(j=0;j<=500;j++) { B[i][j]=A[i][j]; } } lams=fanmifa(); printf("lams=%13.11le\n",lams);//反幂法求模最小的特征值 for(t=1;t<=39;t++)//循环求与miu接近的特征值 { miu=min+t*(max-min)/40; for(i=0;i<=4;i++)//赋值B=(A-miu*I) { for(j=0;j<=500;j++) { if(i==2) B[i][j]=A[i][j]-miu; else B[i][j]=A[i][j]; } } lam=fanmifa()+miu;//用带原点的反幂法求与miu接近的特征值printf("lami%d=%13.11le\n",t,lam); } cond2=fabs(lams/lamk);//计算A的条件数 printf("condA=%13.11le\n",cond2); for(i=0;i<=4;i++)//赋值B=A { for(j=0;j<=500;j++) { B[i][j]=A[i][j]; } } fenjie();//对B进行LU分解，其第三行的元素即为矩阵U的对角线元素，即乘积为A的行列式的值 for(i=0;i<=500;i++) { det=det*B[2][i]; } printf("det(A)=%13.11le\n",det);//计算输出A的行列式 } /********幂法程序********/ double mifa()//幂法的函数 { inti,j,p,s; double u[501],y[501]={0}; double ita,a,beta1=0,beta2=0,b,c,lam; for(i=0;i<=500;i++) u[i]=1; while(1)//开始迭代 { a=0; for(i=0;i<=500;i++) { a=a+u[i]*u[i]; } ita=sqrt(a);//计算2范数 for(i=0;i<=500;i++)//计算向量y(k-1) { y[i]=u[i]/ita; } for(i=0;i<=500;i++) { u[i]=0; } for(i=2;i<=502;i++)//计算u(k) { if(i<4) p=i; else p=4; if(i<500) s=i; else s=500; for(j=i-p;j<=s;j++) { u[i-2]=u[i-2]+B[i-j][j]*y[j]; } } for(i=0;i<=500;i++) beta2=beta2+y[i]*u[i];//计算beta(T) b=fabs(beta2-beta1); c=fabs(beta2); if((b/c)<=sigma)//如果达到精度则结束并返回值 { lam=beta2; goto finish; } else { beta1=beta2; beta2=0; } } finish:return(lam); } /********反幂法程序********/ double fanmifa()//反幂法函数 { inti,k,t,p; double u[501],y[501]={0},d[501]={0}; double ita,a,beta1=0,beta2=0,b,c,lam; for(i=0;i<=500;i++) u[i]=1; fenjie();//对系数矩阵进行LU分解 /*for(i=0;i<=500;i++) printf("%12le ",B[0][i]);*/ for(k=1;;k++)//开始迭代 { a=0; for(i=0;i<=500;i++) { a=a+u[i]*u[i]; } ita=sqrt(a);//计算2范数 for(i=0;i<=500;i++)//计算向量y（k-1） { y[i]=u[i]/ita; } for(i=0;i<=500;i++) { u[i]=0; } d[0]=y[0];//用LU分解解拟三角矩阵的方法求解u(k) for(i=1;i<=500;i++) { t=1; if(i-1>t) t=i-1; else t=t; double temp1=0; for(;t<=i;t++) { temp1+=B[i-t+3][t-1]*d[t-1]; } d[i]=y[i]-temp1; } u[500]=d[500]/B[2][500]; for(i=499;i>=0;i--) { if(i+2>500) p=500; else p=i+2; double temp2=0; for(t=i+1;t<=p;t++) { temp2+=B[i-t+2][t]*u[t]; } u[i]=(d[i]-temp2)/B[2][i]; } for(i=0;i<=500;i++) beta2=beta2+y[i]*u[i]; b=fabs(beta2-beta1); c=fabs(beta2); if((b/c)<=sigma)//如果达到精度则结束并返回值 { lam=1/beta2; goto finish; } else { beta1=beta2; beta2=0; } } finish:return(lam); } /********矩阵分解程序********/ void fenjie()//LU分解函数 { intk,t,j,i,p; for(i=3;i<=4;i++)//对于k=1时对应的矩阵元素先赋值 { B[i][0]=B[i][0]/B[2][0]; } for(k=2;k<=501;k++)//从k=2开始将U中的元素存到系数矩阵中{ if(k+2<501) p=k+2; else p=501; for(j=k;j<=p;j++) { t=1; if(k-2>t) t=k-2; else t=t; if(j-2>t) t=j-2; else t=t; double temp3=0; for(;t<=k-1;t++) { temp3+=B[k-t+2][t-1]*B[t-j+2][j-1]; } B[k-j+2][j-1]=B[k-j+2][j-1]-temp3; } if(k<501) { if(k+2<501) p=k+2; else p=501; for(i=k+1;i<=p;i++) { t=1; if(i-2>t) t=i-2; else t=t; if(k-2>t) t=k-2; else t=t; double temp4=0; for(;t<=k-1;t++) { temp4+=B[i-t+2][t-1]*B[t-k+2][k-1]; } B[i-k+2][k-1]=B[i-k+2][k-1]-temp4; B[i-k+2][k-1]=B[i-k+2][k-1]/B[2][k-1]; } } } } 三、计算结果 计算结果如下图所示 四、编程中遇到的问题 在幂法与反幂法的求解过程中初始值选取不同得到的结果也不同，比如取u0=[1,0,…0]和取u0=[1,1,….1]得到的结果不同，书上说可以任意取初始向量，但这里却不可以，原因未知。 《数值分析B》大作业一 SY1103120 朱舜杰 一．算法设计方案： 1.矩阵A的存储与检索 将带状线性矩阵A[501][501]转存为一个矩阵MatrixC[5][501] . 由于C语言中数组角标都是从0开始的，所以在数组MatrixC[5][501]中检索A的带内元素a ij的方法是： A的带内元素a ij=C中的元素c i-j+2,j 2.求解λ1，λ501，λs ①首先分别使用幂法和反幂法迭代求出矩阵按摸最大和最小的特征值λmax和λmin。λmin即为λs； 如果λmax>0,则λ501=λmax；如果λmax<0,则λ1=λmax。 ②使用带原点平移的幂法（mifa（）函数），令平移量p=λmax，求 出对应的按摸最大的特征值λ，max， 如果λmax>0,则λ1=λ，max+p；如果λmax<0,则λ501=λ，max+p。 3.求解A的与数μk=λ1+k（λ501-λ1）/40的最接近的特征值λik （k=1,2，…，39）。 使用带原点平移的反幂法，令平移量p=μk，即可求出与μk最接近的特征值λik。 4.求解A的（谱范数）条件数cond（A）2和行列式d etA。 ①cond（A）2=|λ1/λn|，其中λ1和λn分别是矩阵A的模最大和 最小特征值。 ②矩阵A的行列式可先对矩阵A进行LU分解后，detA等于U所有对角线上元素的乘积。 二．源程序 #include 数值分析计算实习题 第二章 2-1 程序： clear;clc; x1=[0.2 0.4 0.6 0.8 1.0]; y1=[0.98 0.92 0.81 0.64 0.38]; n=length(y1); c=y1(:); for j=2:n %求差商 for i=n:-1:j c(i)=(c(i)-c(i-1))/(x1(i)-x1(i-j+1)); end end syms x df d; df(1)=1;d(1)=y1(1); for i=2:n %求牛顿差值多项式 df(i)=df(i-1)*(x-x1(i-1)); d(i)=c(i)*df(i); end disp('4次牛顿插值多项式'); P4=vpa(collect((sum(d))),5) %P4即为4次牛顿插值多项式,并保留小数点后5位数pp=csape(x1,y1, 'variational');%调用三次样条函数 q=pp.coefs; disp('三次样条函数'); for i=1:4 S=q(i,:)*[(x-x1(i))^3;(x-x1(i))^2;(x-x1(i));1]; S=vpa(collect(S),5) end x2=0.2:0.08:1.08; dot=[1 2 11 12]; figure ezplot(P4,[0.2,1.08]); hold on y2=fnval(pp,x2); x=x2(dot); y3=eval(P4); y4=fnval(pp,x2(dot)); plot(x2,y2,'r',x2(dot),y3,'b*',x2(dot),y4,'co'); title('4次牛顿插值及三次样条'); 结果如下： 4次牛顿插值多项式 P4 = - 0.52083*x^4 + 0.83333*x^3 - 1.1042*x^2 + 0.19167*x + 0.98 三次样条函数 《数值分析》计算实习题目 第一题： 1. 算法设计方案 （1）1λ，501λ和s λ的值。 1)首先通过幂法求出按模最大的特征值λt1，然后根据λt1进行原点平移求出另一特征值λt2，比较两值大小，数值小的为所求最小特征值λ1，数值大的为是所求最大特征值λ501。 2)使用反幂法求λs ，其中需要解线性方程组。因为A 为带状线性方程组，此处采用LU 分解法解带状方程组。 （2）与140k λλμλ-5011=+k 最接近的特征值λik 。 通过带有原点平移的反幂法求出与数k μ最接近的特征值 λik 。 （3）2cond(A)和det A 。 1）1=n λλ2cond(A)，其中1λ和n λ分别是按模最大和最小特征值。 2）利用步骤（1）中分解矩阵A 得出的LU 矩阵，L 为单位下三角阵,U 为上三角阵，其中U 矩阵的主对角线元素之积即为det A 。 由于A 的元素零元素较多，为节省储存量，将A 的元素存为6×501的数组中，程序中采用get_an_element()函数来从小数组中取出A 中的元素。 2.全部源程序 #include 插值法 1.下列数据点的插值 x 0 1 4 9 16 25 36 49 64 y 0 1 2 3 4 5 6 7 8 可以得到平方根函数的近似，在区间[0,64]上作图. (1)用这9个点作8次多项式插值Ls(x). (2)用三次样条(第一边界条件)程序求S(x). 从得到结果看在[0,64]上，哪个插值更精确；在区间[0,1]上，两种插值哪个更精确？ 解：(1)拉格朗日插值多项式，求解程序如下 syms x l; x1=[0 1 4 9 16 25 36 49 64]; y1=[0 1 2 3 4 5 6 7 8]; n=length(x1); Ls=sym(0); for i=1:n l=sym(y1(i)); for k=1:i-1 l=l*(x-x1(k))/(x1(i)-x1(k)); end for k=i+1:n l=l*(x-x1(k))/(x1(i)-x1(k)); end Ls=Ls+l; end Ls=simplify(Ls) %为所求插值多项式Ls(x). 输出结果为 Ls = -24221063/63504000*x^2+95549/72072*x-1/3048192000*x^8-2168879/435456000 *x^4+19/283046400*x^7+657859/10886400*x^3+33983/152409600*x^5-13003/2395008 000*x^6 (2)三次样条插值，程序如下 x1=[0 1 4 9 16 25 36 49 64]; y1=[0 1 2 3 4 5 6 7 8]; x2=[0:1:64]; y3=spline(x1,y1,x2); p=polyfit(x2,y3,3); %得到三次样条拟合函数 S=p(1)+p(2)*x+p(3)*x^2+p(4)*x^3 %得到S(x) 输出结果为： S = 23491/304472833/8*x+76713/*x^2+6867/42624*x^3 (3)在区间[0,64]上，分别对这两种插值和标准函数作图， plot(x2,sqrt(x2),'b',x2,y2,'r',x2,y3,'y') 蓝色曲线为y=函数曲线,红色曲线为拉格朗日插值函数曲线，黄色曲线为三次样条插值曲线 010203040506070 -200 20 40 60 80100 可以看到蓝色曲线与黄色曲线几乎重合，因此在区间[0,64]上三次样条插值更精确。 在[0,1]区间上由上图看不出差别，不妨代入几组数据进行比较 ，取x4=[0:0.2:1] 数值分析第三次大作业 一、算法的设计方案： （一）、总体方案设计： x y当作已知量代入题目给定的非线性方程组，求（1）解非线性方程组。将给定的(,) i i 得与(,)i i x y 相对应的数组t[i][j],u[i][j]。 （2）分片二次代数插值。通过分片二次代数插值运算，得到与数组t[11][21],u[11][21]]对应的数组z[11][21]，得到二元函数z=(,)i i f x y 。 （3）曲面拟合。利用x[i],y[j],z[11][21]建立二维函数表，再根据精度的要求选择适当k 值，并得到曲面拟合的系数矩阵C[r][s]。 （4）观察和(,)i i p x y 的逼近效果。观察逼近效果只需要重复上面（1）和（2）的过程，得到与新的插值节点(,)i i x y 对应的(,)i i f x y ，再与对应的(,)i i p x y 比较即可，这里求解 (,)i i p x y 可以直接使用（3）中的C[r][s]和k 。 （二）具体算法设计： （1）解非线性方程组 牛顿法解方程组()0F x =的解* x ，可采用如下算法： 1）在* x 附近选取(0) x D ∈，给定精度水平0ε>和最大迭代次数M 。 2）对于0,1, k M =执行 ① 计算() ()k F x 和()()k F x '。 ② 求解关于() k x ?的线性方程组 () ()()()()k k k F x x F x '?=- ③ 若() () k k x x ε∞∞ ?≤，则取*()k x x ≈，并停止计算；否则转④。 ④ 计算(1) ()()k k k x x x +=+?。 ⑤ 若k M <，则继续，否则，输出M 次迭代不成功的信息，并停止计算。 （2）分片双二次插值 给定已知数表以及需要插值的节点，进行分片二次插值的算法： 设已知数表中的点为： 00(0,1,,) (0,1,,)i j x x ih i n y y j j m τ=+=???=+=?? ，需要插值的节点为(,)x y 。 1) 根据(,)x y 选择插值节点(,)i j x y ： 若12h x x ≤+ 或12 n h x x ->-，插值节点对应取1i =或1i n =-， 目标：使用带双步位移的QR 分解法求矩阵10*10[]ij A a =的全部特征值，并对其中的每一个实特征值求相应的特征向量。已知：sin(0.50.2)() 1.5cos( 1.2)(){i j i j ij i j i j a +≠+== (i,j=1,2, (10) 算法： 以上是程序运作的逻辑，其中具体的函数的算法，大部分都是数值分析课本上的逻辑，在这里特别写出矩阵A 的实特征值对应的一个特征向量的求法： ()[]()() []()[]()111111I 00000 i n n n B A I gause i n Q A I u Bu u λλ-?-?-=-?-?? ?-=????→=??????→= ?? ? 选主元的消元 检查知无重特征值 由于=0i A I λ- ，因此在经过选主元的高斯消元以后，i A I λ- 即B 的最后一行必然为零，左上方变 为n-1阶单位矩阵[]()()11I n n -?-，右上方变为n-1阶向量[]()11n Q ?-，然后令n u 1=-，则 ()1,2,,1j j u Q j n ==???-。 这样即求出所有A所有实特征值对应的一个特征向量。 #include 第四章： 1、（1）：复合梯形 建立m文件： function t=natrapz(fname,a,b,n) h=(b-a)/n; fa=feval(fname,a);fb=feval(fname,b); f=feval(fname,a+h:h:b-h+0.001*h); t=h*(0.5*(fa+fb)+sum(f)); 输入： >> syms x >> f=inline('sqrt(x).*log(x);'); >> natrapz(f,eps,1,10) 输出： ans = -0.417062831779470 输入： >> syms x >> f=inline('sqrt(x).*log(x);'); >> natrapz(f,eps,1,100) 输出： ans = -0.443117908008157 输入： >> syms x >> f=inline('sqrt(x).*log(x);'); >> natrapz(f,eps,1,1000) 输出： ans = -0.444387538997162 复合辛普森 建立m文件： function t=comsimpson(fname,a,b,n) h=(b-a)/n; fa=feval(fname,a);fb=feval(fname,b); f1=feval(fname,a+h:h:b-h+0.001*h); f2=feval(fname,a+h/2:h:b-h+0.001*h); t=h/6*(fa+fb+2*sum(f1)+4*sum(f2)); 输入： >> syms x >> f=inline('sqrt(x).*log(x);'); >> format long； >>comsimpson(f,eps,1,10) 输出： ans = -0.435297890074689 输入： >>syms x >>f=inline('sqrt(x).*log(x);'); >>comsimpson(f,eps,1,100) 输出： ans = -0.444161178415673 输入： >>syms x >>f=inline('sqrt(x).*log(x);'); >>comsimpson(f,eps,1,1000) 输出： ans = -0.444434117614180 （2）龙贝格 建立m文件： function [RT,R,wugu,h]=Romberg(fun,a,b,wucha,m) %RT是龙贝格积分表 %R是数值积分值 %wugu是误差估计 %h是最小步长 %fun是被积函数 %a b是积分下、上限 数值分析第二次大作业 史立峰 SY1505327 一、 方案 （1）利用循环结构将sin(0.50.2)() 1.5cos( 1.2)() {i j i j ij i j i j a +≠+==(i,j=1,2,……,10)进行赋值，得到需要变换的 矩阵A ； （2）然后，对矩阵A 利用Householder 矩阵进行相似变换，把A 化为上三角矩阵A (n-1)。 对A 拟上三角化，得到拟上三角矩阵A (n-1)，具体算法如下： 记A(1)=A ，并记A(r)的第r 列至第n 列的元素为()n r r j n i a r ij ,,1,;,,2,1) ( +==。 对于2,,2,1-=n r 执行 1. 若 ()n r r i a r ir ,,3,2) ( ++=全为零，则令A(r+1) =A(r),转5；否则转2。 2. 计算 () ∑+== n r i r ir r a d 1 2 )( ()( )r r r r r r r r r r d c a d a c ==-=++则取,0sgn ) (,1)(,1若 )(,12r r r r r r a c c h +-= 3. 令 () n T r nr r r r r r r r r R a a c a u ∈-=++) ()(,2)(,1,,,,0,,0 。 4. 计算 r r T r r h u A p /)(= r r r r h u A q /)(= r r T r r h u p t /= r r r r u t q -=ω T r r T r r r r p u u A A --=+ω)()1( 5. 继续。 （3）使用带双步位移的QR 方法计算矩阵A (n-1)的全部特征值，也是A 的全部特征值，具体算法如下： 1. 给定精度水平0>ε和迭代最大次数L 。 2. 记n n ij n a A A ?-==][) 1()1()1(，令n m k ==,1。 1． 设?+=1 05dx x x I n n ， （1） 由递推公式n I I n n 1 51+-=-，从0I 的几个近似值出发，计算20I ； 解：易得：0I =ln6-ln5=0.1823, 程序为： I=0.182; for n=1:20 I=(-5)*I+1/n; end I 输出结果为：20I = -3.0666e+010 （2） 粗糙估计20I ，用n I I n n 51 5111+- =--，计算0I ； 因为 0095.05 6 0079.01020 201 020 ≈<<≈??dx x I dx x 所以取0087.0)0095.00079.0(2 1 20=+= I 程序为：I=0.0087; for n=1:20 I=(-1/5)*I+1/(5*n); end I 0I = 0.0083 （3） 分析结果的可靠性及产生此现象的原因(重点分析原因)。 首先分析两种递推式的误差；设第一递推式中开始时的误差为000I I E '-=，递推过程的舍入误差不计。并记n n n I I E '-=，则有01)5(5E E E n n n -==-=-Λ。因为=20E 20020)5(I E >>-，所此递推式不可靠。而在第二种递推式中n n E E E )5 1(5110-==-=Λ，误差在缩小， 所以此递推式是可靠的。出现以上运行结果的主要原因是在构造递推式过程中，考虑误差是否得到控制， 即算法是否数值稳定。 2． 求方程0210=-+x e x 的近似根，要求4 1105-+?<-k k x x ，并比较计算量。 （1） 在[0，1]上用二分法； 程序：a=0;b=1.0; while abs(b-a)>5*1e-4 c=(b+a)/2; 北京航空航天大学 数值分析大作业八 学院名称自动化 专业方向控制工程 学号 学生姓名许阳 教师孙玉泉 日期2014 年11月26 日 一．题目 关于x , y , t , u , v , w 的方程组(A.3) ???? ?? ?=-+++=-+++=-+++=-+++79 .0sin 5.074.3cos 5.007.1cos sin 5.067.2cos 5.0y w v u t x w v u t y w v u t x w v u t (A.3) 以及关于z , t , u 的二维数表（见表A-1）确定了一个二元函数z =f (x , y )。 表A-1 二维数表 t z u 0 0.4 0.8 1.2 1.6 2 0 -0.5 -0.34 0.14 0.94 2.06 3.5 0.2 -0.42 -0.5 -0.26 0.3 1.18 2.38 0.4 -0.18 -0.5 -0.5 -0.18 0.46 1.42 0.6 0.22 -0.34 -0.58 -0.5 -0.1 0.62 0.8 0.78 -0.02 -0.5 -0.66 -0.5 -0.02 1.0 1.5 0.46 -0.26 -0.66 -0.74 -0.5 1. 试用数值方法求出f (x , y ) 在区域}5.15.0,8.00|), {≤≤≤≤=y x y x D （上的近似表达式 ∑∑===k i k j s r rs y x c y x p 00 ),( 要求p (x , y )以最小的k 值达到以下的精度 ∑∑==-≤-=10020 7210)],(),([i j i i i i y x p y x f σ 其中j y i x i i 05.05.0,08.0+==。 2. 计算),(),,(* ***j i j i y x p y x f (i =1,2,…,8 ; j =1,2,…,5) 的值，以观察p (x , y ) 逼 近f (x , y )的效果，其中j y i x j i 2.05.0,1.0**+==。 《数值分析》计算实习题 姓名： 学号： 班级： 第二章 1、程序代码 Clear;clc; x1=[0.2 0.4 0.6 0.8 1.0]; y1=[0.98 0.92 0.81 0.64 0.38]; n=length(y1); c=y1(:); for j=2:n %求差商 for i=n:-1:j c(i)=(c(i)-c(i-1))/(x1(i)-x1(i-j+1)); end end syms x df d; df(1)=1;d(1)=y1(1); for i=2:n %求牛顿差值多项式 df(i)=df(i-1)*(x-x1(i-1)); d(i)=c(i-1)*df(i); end P4=vpa(sum(d),5) %P4即为4次牛顿插值多项式,并保留小数点后5位数 pp=csape(x1,y1, 'variational');%调用三次样条函数 q=pp.coefs; q1=q(1,:)*[(x-.2)^3;(x-.2)^2;(x-.2);1]; q1=vpa(collect(q1),5) q2=q(1,:)*[(x-.4)^3;(x-.4)^2;(x-.4);1]; q2=vpa(collect(q2),5) q3=q(1,:)*[(x-.6)^3;(x-.6)^2;(x-.6);1]; q3=vpa(collect(q3),5) q4=q(1,:)*[(x-.8)^3;(x-.8)^2;(x-.8);1]; q4=vpa(collect(q4),5)%求解并化简多项式 2、运行结果 P4 = 0.98*x - 0.3*(x - 0.2)*(x - 0.4) - 0.625*(x - 0.2)*(x - 0.4)*(x - 0.6) - 0.20833*(x - 0.2)*(x - 0.4)*(x - 0.8)*(x - 0.6) + 0.784 q1 = - 1.3393*x^3 + 0.80357*x^2 - 0.40714*x + 1.04 数值分析——数值积分实习题 管理科学与工程学院 学号：1120140500 姓名：彭洋洋 一、计算实习题 1. 用不同数值方法计算积分：04 9 xdx =-?. （1）取不同的步长h ，分别用复合梯形及复合辛普森求积计算积分，给出误差中关于h 的函数，并与积分精确值比较两个公式的精度，是否存在一个最小的h ，使得精度不能再被改善？ （2）用龙贝格求积计算完成问题（1） （3）用自适应辛普森积分，使其精度达到10-4 解答：（1）取不同的步长，采用不同的公式，比较精度过程如下： 1.1 复合梯形公式及复合辛普森公式求解 复合梯形公式：1 1 *[()2()()]2n n k k h T f a f x f b -==++∑ 误差关于h 的函数：2 (2)()**()12 n a b R f h f ξ-= 复合辛普森公式：11 1/201 *[()4()2()()]6n n n k k k k h S f a f x f x f b --+===+++∑∑ 误差关于h 的函数：4(4)()*(/2)*()180 n a b R f h f η-= 1.2 复合梯形公式及复合辛普森公式Matlab 程序 （2）用龙贝格求积计算完成问题（1） 2.1 龙贝格求积算法 龙贝格求积公式也称为逐次分半加速法。它是在梯形公式、辛普森公式和柯特斯公式之间的关系的基础上，构造出一种加速计算积分的方法。作为一种外推算法，它在不增加计算量的前提下提高了误差的精度。 24133n n n S T T =- 21611515n n n C S S =- 26416363 n n n R C C =- 1221/20 1()22n n n k k h T T f x -+==+∑ ()(1)()11(4*)/(41)k m k k m m m m T T T +--=-- 1,2,.. k = 2.2 龙贝格求积Matlab 程序 数值分析计算实习题（二） SY1004114 全昌彪一：算法设计方案概述： 本题采用fortran90语言编写程序，依据题目要求，采用带双步位移QR分解法求出所给矩阵的所有特征值，并求出相应于其实特征值的特征向量，以及相关需要给出的中间结果。 1、矩阵的A的初始化（赋值）：利用子函数initial（a,n）来实现，返回n×n 维二维数组a。 2、A矩阵的拟上三角化：利用子函数hessenberg(a,n)，在对矩阵进行QR分解前进行拟上三角化，这样可以提高计算效率，减少计算量，返回A矩阵的相似矩阵Hessenberg阵A(n-1)。 3、对A(n-1)进行带双步位移QR分解得出Ｃｍ及A矩阵的所有特征值，这一步利用了两个子函数eigenvalue(a,n,lamda,lamdai)和qrresolve(b,c,m) 带双步位移QR分解可以加速收敛。每次QR分解前先进行判断，若可以直接得到矩阵的特征值，则对矩阵直接降阶处理；若不可以，则进行QR分解，这样就进一步减少了计算量，提高了计算效率。考虑到矩阵A可能有复特征值，采用两个一维数组lamda(n)及lamdai(n)分别存储其实部和虚部。在双步位移处理及降阶过程中，被分解的矩阵Ａk(m ×m)及中间矩阵M k(m×m)的维数随m不断减少而降阶，于是引入了动态矩阵C(m×m)和B(m×m)分别存储，在使用前，先声明分配内存，使用结束后立即释放内存。返回A(n-1)经双步位移QR分解后的矩阵及Ａ矩阵的所有特征值。 4、特征向量的求解：采用子函数eigenvector(a,lamda)实现求解A矩阵的属于实特征值的特征向量。核心算法为高斯列主元消去法，(A-λI)x=b,b=0,回代过程令x(10)=1，即可求出对应于每一实特征值的特征向量的各个元素。 5、相关输出结果：所有数据均采用e型输出，数据保留到小数点后12位。对于Ａ矩阵的拟上三角化结果集双步位移QR分解结果比较庞大，为了更好的显示，还采用了ｆ型输出，保留５位小数。所有计算精度水平Ｅ＝10-12。 二：算法fortran源程序 !此函数用于给A赋值，即初始化A矩阵 subroutine initial(a,n) integer :: i,j dimension a(n,n) double precision a do i=1,n do j=1,n if (i==j) then a(i,j)=1.5*cos(i+1.2*j) else a(i,j)=sin(0.5*i+0.2*j) end if end do end do 《数值分析A》计算实习题目第一题 一．算法设计方案： 1.矩阵A的存储与检索 将带状线性矩阵A[501][501]转存为一个矩阵MatrixC[5][501] . 由于C语言中数组角标都是从0开始的，所以在数组MatrixC[5][501]中检索A的带内元素a ij的方法是： A的带内元素a ij=C中的元素c i-j+2,j 2.求解λ1，λ501，λs ①首先分别使用幂法和反幂法迭代求出矩阵按摸最大和最小的特征值λmax和λmin。λmin即为λs； 如果λmax>0,则λ501=λmax；如果λmax<0,则λ1=λmax。 ②使用带原点平移的幂法（mifa（）函数），令平移量p=λmax，求出对应的按摸最大的特征值λ，max， 如果λmax>0,则λ1=λ，max+p；如果λmax<0,则λ501=λ，max+p。 3.求解A的与数μk=λ1+k（λ501-λ1）/40的最接近的特征值λik （k=1,2，…，39）。 使用带原点平移的反幂法，令平移量p=μk，即可求出与μk最接近的特征值λik。 4.求解A的（谱范数）条件数cond（A）2和行列式d etA。 ①cond（A）2=|λ1/λn|，其中λ1和λn分别是矩阵A的模最大和最小特征值。 ②矩阵A的行列式可先对矩阵A进行LU分解后，detA等于U所有 对角线上元素的乘积。 二．源程序(VS2010环境下，C++语言) #include 数值分析计算实习题一 学号： ： 院系： 2015年11月5日 一、分析 1.1算法分析 题目要求求出： 1）特征值从小到大排列的最小特征值1λ和最大特征值501λ。 2）特征值中模最小的特征值s λ。 3）靠近一组数k μ的一组特征值k i λ。 4）矩阵A 的条件数cond(A)2。 5）行列式detA 。 解决方法： 1）若将所有行列式按模的大小排列则模最大的特征值一定是1λ和501λ中的一个，因此利用幂法求出模最大的特征值1m λ。然后利用带原点平移的幂法，将系数矩阵变为1m A I λ-即将所有特征值都减去1m λ，则特征值按大小顺序排列的次序不变，模最大的特征值依然在整个排列的两端，再用一次幂法得到模最大的特征值 21=m m λλλ-，其中λ为带原点平移的幂法求出的特征值，最后两个特征值1m λ、2m λ比较大小，大的为501λ，小的为1λ。 2）因为s λ为按模最小的特征值，因此用反幂法可求的其特征值。 3）因为k i λ靠近数k μ，因此k i k λμ-一定是所有的k λμ-中模最小的，因此可利用带原点平移的反幂法求出特征值k i λ，此时的系数矩阵变为k A I μ-。 4）条件数cond(A)2为模最小的特征值与模最大的特征值的比的绝对值，因此利用1和2中求出的1m λ和s λ可解出条件数。 5）可对矩阵A 进行LU 分解，即A LU =则det()det()det()A L U =?，又因为矩阵L 对角线元素为1，则det()L =1，所以det()det()A U =，U 为上三角阵，行列式为对角线元素的乘积，因此可得A 的行列式。 1.2程序分析 一、算法设计方案 1.使用牛顿迭代法，对原题中给出的i x i 08.0=，j y j 05.05.0+=， (010 ,020i j ≤≤≤≤)的11*21组j i y x ,分别求出原题中方程组的一组解，于是得到一组和i i y x ,对应的j i t u ,。 2.对于已求出的j i t u ,，使用分片二次代数插值法对原题中关于u t z ,,的数表进行插值得到 ij z 。于是产生了z=f(x,y)的11*21个数值解。 3.从k=1开始逐渐增大k 的值，并使用最小二乘法曲面拟合法对z=f(x,y)进行拟合，得到每次的σ,k 。当7 10-<σ时结束计算，输出拟合结果。 4.计算)5,,2,1,8,,2,1)(,(),,(* ***???=???=j i y x p y x f j i j i 的值并输出结果，以观察),(y x p 逼近),(y x f 的效果。其中j y i x j i 2.05.0,1.0* *+==。 二、算法实现方案 1、求(,)f x y ： （1）Newton 法解非线性方程组 0.5cos 2.670.5sin 1.07(1)0.5cos 3.740.5sin 0.79 t u v w x t u v w y t u v w x t u v w y +++-=??+++-=? ? +++-=??+++-=?， 其中，t, u, v ,w 为待求的未知量，x, y 为代入的已知量。 设(,,,)T t u v w ξ=，给定精度水平12110ε-=和最大迭代次数M ，则解该线性方程组的迭代格式为： *(0)(0)(0)(0)(0)(k+1) ()()1()(,,,)()()0,1,T k k k t u v w F F k ξξξ ξξξ-?=?'=-??= ? 在附近选取初值， 迭代终止条件为()(1) () 1/k k k ξξ ξε-∞ ∞ -≤，若k M >时仍未达到迭代精度，则迭代计算失 败。 其中，雅可比矩阵 0.5*cos(t) + u + v + w - x - 2.67t + 0.5*sin(u) + v + w - y - 1.07()0.5*t + u + cos(v) + w - x - 3.74t + 0.5*u + v + sin(w) - y - 0.79F ξ???? ? ?=?????? ， 航空航天大学 《数值分析》计算实习报告 第一大题 学院：自动化科学与电气工程学院 专业：控制科学与工程 学生姓名： 学号： 教师： 电话： 完成日期： 2015年11月6日 航空航天大学 Beijing University of Aeronautics and Astronautics 实习题目： 第一题 设有501501?的实对称矩阵A ， ??? ???? ?????????=5011A a b c b c c b c b a 其中，064.0,16.0),501,,2,1(64.0)2.0sin()024.064.1(1 .0-==???=--=c b i e i i a i i 。矩阵A 的特征值为)501,,2,1(???=i i λ，并且有 ||min ||,501 150121i i s λλλλλ≤≤=≤???≤≤ 1.求1λ，501λ和s λ的值。 2.求A 的与数40 1 5011λλλμ-+=k k 最接近的特征值)39,,2,1(???=k k i λ。 3.求A 的(谱数)条件数2)A (cond 和行列式detA 。 说明： 1.在所用的算法中，凡是要给出精度水平ε的，都取12-10=ε。 2.选择算法时，应使矩阵A 的所有零元素都不储存。 3.打印以下容： （1）全部源程序； （2）特征值），，39,...,2,1(,s 5011=k k i λλλλ以及A det ,)A (cond 2的值。 4.采用e 型输出实型数，并且至少显示12位有效数字。 一、算法设计方案 1、求1λ，501λ和s λ的值。 由于||min ||,501 150121i i s λλλλλ≤≤=≤???≤≤，可知绝对值最大特征值必为1λ和501 λ其中之一，故可用幂法求出绝对值最大的特征值λ，如果λ=0，则1λ=λ，否则 501λ=λ。将矩阵A 进行一下平移： I -A A'λ= （1） 对'A 用幂法求出其绝对值最大的特征值'λ，则A 的另一端点特征值1λ或501λ为'λ+λ。 s λ为按模最小特征值，||min ||501 1i i s λλ≤≤=，可对A 使用反幂法求得。 2、求A 的与数40 1 5011λλλμ-+=k k 最接近的特征值)39,...,2,1(=k k i λ。 计算1)1,2,...,50=(i i λ-k μ，其模值最小的值对应的特征值k λ与k μ最接近。因此对A 进行平移变换： )39,,2,1k -A A k k ==（I μ （2） 对k A 用反幂法求得其模最小的特征值'k λ，则k λ='k λ+k μ。 3、求A 的(谱数)条件数2)(A cond 和行列式detA 。 由矩阵A 为非奇异对称矩阵可得： | | )(min max 2λλ=A cond （3） 其中max λ为按模最大特征值，min λ为按模最小特征值，通过第一问我们求得的λ和s λ可以很容易求得A 的条件数。 在进行反幂法求解时，要对A 进行LU 分解得到。因L 为单位下三角阵，行 列式为1，U 为上三角阵，行列式为主对角线乘积，所以A 的行列式等于U 的行列式，为U 的主对角线的乘积。 北京航空航天大学2020届研究生 《数值分析》实验作业 第九题 院系：xx学院 学号： 姓名： 2020年11月 Q9:方程组A.4 一、 算法设计方案 （一）总体思路 1.题目要求∑∑=== k i k j s r rs y x c y x p 00 ),(对f(x, y) 进行拟合，可选用乘积型最小二乘拟合。 ),(i i y x 与),(i i y x f 的数表由方程组与表A-1得到。 2.),(* * j i y x f 与1使用相同方法求得，),(* * j i y x p 由计算得出的p(x,y)直接带入),(* * j i y x 求得。 1. ),(i i y x 与),(i i y x f 的数表的获得 对区域D ={ （x,y)|1≤x ≤1.24,1.0≤y ≤1.16}上的f (x , y )值可通过xi=1+0.008i ，yj=1+0.008j ，得到),(i i y x 共31×21组。将每组带入A4方程组，即可获得五个二元函数组，通过简单牛顿迭代法求解这五个二元数组可获得z1~z5有关x,y 的表达式。再将 ),(i i y x 分别带入z1~z5表达式即可获得f(x,y)值。 2.乘积型最小二乘曲面拟合 2.1使用乘积型最小二乘拟合，根据k 值不用，有基函数矩阵如下： ????? ??=k i i k x x x x B 0000 ， ????? ??=k j j k y y y y G 0000 数表矩阵如下： ???? ? ? ?=),(),(),(),(0000j i i j y x f y x f y x f y x f U 记C=[rs c ]，则系数rs c 的表达式矩阵为： 11-)(-=G G UG B B B C T T T ）（ 通过求解如下线性方程，即可得到系数矩阵C 。 UG B G G C B B T T T =)()( 2.2计算),(),,(* ***j i j i y x p y x f (i =1,2,…,31 ; j =1,2,…,21) 的值 ),(**j i y x f 的计算与),(j i y x f 相同。将),(**j i y x 代入原方程组，求解响应) ,(* *ij ij u t 进行分片双二次插值求得),(**j i y x f 。),(* *j i y x p 的计算则可以直接将),(**j i y x 代入所求p(x,y)。 二、 源程序 ********* 第三次数值分析大作业Q9************ integer::i, j, K1, L1, n, m dimension X(31), Y(21), T(6), U(6), Z(6, 6), UX(11, 21), TY(11, 21), FXY(11, 21), C(6, 6) dimension z1(31, 21), z2(31, 21), z3(31, 21), z4(31, 21), z5(31, 21) dimension X1(8), Y1(5), FXY1(8, 5), PXY1(8, 5), UX1(8, 5), TY1(8, 5)北航数值分析大作业一
数值分析第五版计算实习题
北航数值分析大作业第一题幂法与反幂法
数值分析计算实习题
北航数值分析报告第三次大作业
北航数值分析大作业第二题精解
数值分析(第五版)计算实习题第四章作业
北航数值分析大作业第二题
(完整版)数值计算方法上机实习题答案
北航数值分析报告大作业第八题
数值分析计算实习题
数值分析课程实验设计——数值积分实习题
数值分析计算实习题(二)
北航数值分析课程第一次大作业讲解
数值分析计算实习作业一
北航数值分析大作业3
北航数值分析计算实习报告一
BUAA数值分析大作业三
相关主题
文本预览