数值分析计算实习题

《数值分析》计算实习题二算法设计方案1．主要计算步骤：计算函数f(x,y)在拟合所需的节点处的函数值。

将各拟合节点（x i，y j）分别带入非线性方程组0.5 cos t + u + v + w – x = 2.67t + 0.5 sin u + v + w – y = 1.070.5t + u + cos v + w – x =3.74t + 0.5u + v + sin w – y =0.79解非线性方程组得解向量（t ij，u ij，v ij，w ij）。

对数表z(t,u)进行分片二次代数插值，求得对应（t ij，u ij）处的值，即为f(x i，y j) 的值。

对上述拟合节点分别进行x，y最高次数为k（k=0,1,2,3…）次的多项式拟合。

每次拟合后验证误差大小，直到满足要求。

2．求解非线性方程组选择Newton迭代法，迭代过程中需要求解线性方程组，选择选主元的Doolittle分解法。

3．对z(t，u)进行插值选择分片二次插值。

4．拟合基函数φr(x)ψs(y)选择为φr(x)=x r，ψs(y)=y s。

拟合系数矩阵c通过连续两次解线性方程组求得。

一．源程序#include "stdio.h"#include "stdlib.h"#include "math.h"void Doolittle(double *A,int n,int *M)//功能说明：对n阶矩阵A进行选主元的Doolittle分解//参数说明：A：欲进行分解的方阵，同时也是返回参数，分解后的结果// 存储于A中// n：方阵A的维数// M；（返回参数）n维向量，记录选主元过程中行交换的次序{int i,j,k,t;double *s;double Maxs,temp;s=(double*) calloc(n,sizeof(double));for(k=0;k<n;k++){for(i=k;i<n;i++){s[i]=A[i*n+k];for(t=0;t<k;t++) s[i]-= A[i*n+t] * A[t*n+k];}Maxs=abs(s[k]); M[k]=k;for(i=k+1;i<n;i++){if(Maxs<abs(s[i])){Maxs=abs(s[i]);M[k]=i;}}if(M[k]!=k){for(t=0;t<n;t++){temp=A[k*n+t];A[k*n+t]=A[M[k]*n+t];A[M[k]*n+t]=temp;}temp=s[k];s[k]=s[M[k]];s[M[k]]=temp;}if(Maxs<(1e-14)){s[k]=1e-14;printf("%.16e方阵奇异\n",Maxs);}A[k*n+k]=s[k];for(j=k+1;(j<n)&&(k<n-1);j++){for(t=0;t<k;t++) A[k*n+j]-=A[k*n+t]*A[t*n+j];A[j*n+k]=s[j]/A[k*n+k];}}}void Solve_LUEquation(double* A,int n,double* b,double* x) //功能说明：解方程LUx＝b，其中L、U共同存储在A中//参数说明：A：经Doolittle分解后的方阵// n：方阵A的维数// b：方程组的右端向量// x：（返回参数）方程组的解向量{int i,t;for(i=0;i<n;i++){x[i]=b[i];for(t=0;t<i;t++) x[i]-=A[i*n+t]*x[t];}for(i=n-1;i>-1;i--){for(t=i+1;t<n;t++) x[i]-=A[i*n+t]*x[t];x[i]/=A[i*n+i];}}void Transpose(double *A,int m,int n,double* AT)//功能说明：求m×n阶矩阵A的转置AT//参数说明：A：已知m×n阶矩阵// m：A的行数// n：A的列数// AT：（返回参数）A的转置矩阵（n×m）{int i,j;for(i=0;i<m;i++)for(j=0;j<n;j++) AT[j*m+i]=A[i*n+j];}void Solve_LEquation(double* A,int n,double* B,double* x,int m) //功能说明：解线性方程组Ax＝B，该函数可对系数矩阵相同// 而右端向量不同的多个方程组同时求解。

数值分析计算实习题-CAL-FENGHAI-(2020YEAR-YICAI)_JINGBIAN《数值分析》计算实习题姓名：学号：班级：第二章1、程序代码Clear;clc;x1=[ ];y1=[ ];n=length(y1);c=y1(:);for j=2:n %求差商for i=n:-1:jc(i)=(c(i)-c(i-1))/(x1(i)-x1(i-j+1));endendsyms 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);endP4=vpa(sum(d),5) %P4即为4次牛顿插值多项式,并保留小数点后5位数pp=csape(x1,y1, 'variational');%调用三次样条函数q=;q1=q(1,:)*[^3;^2;;1];q1=vpa(collect(q1),5)q2=q(1,:)*[^3;^2;;1];q2=vpa(collect(q2),5)q3=q(1,:)*[^3;^2;;1];q3=vpa(collect(q3),5)q4=q(1,:)*[^3;^2;;1];q4=vpa(collect(q4),5)%求解并化简多项式2、运行结果P4 =*x - *(x - *(x - - *(x - *(x - *(x - - *(x - *(x - *(x - *(x - + q1 =- *x^3 + *x^2 - *x +q2 =- *x^3 + *x^2 - *x + q3 =- *x^3 + *x^2 - *x + q4 =- *x^3 + *x^2 - *x +3、问题结果4次牛顿差值多项式4()P x = *x - *(x - *(x - - *(x - *(x - *(x - - *(x - *(x - *(x - *(x - +三次样条差值多项式()Q x0.10.20.30.40.50.60.70.80.910.40.50.60.70.80.911.1323232321.33930.803570.40714 1.04,[0.2,0.4]1.3393 1.60710.88929 1.1643,[0.4,0.6]1.3393 2.4107 1.6929 1.4171,[0.6,0.8]1.3393 3.21432.8179 1.8629,[0.8,1.0]x x x x x x x x x x x x x x x x ⎧-+-+∈⎪-+-+∈⎪⎨-+-+∈⎪⎪-+-+∈⎩第三章1、程序代码Clear;clc; x=[0 1]; y=[1 ];p1=polyfit(x,y,3)%三次多项式拟合 p2=polyfit(x,y,4)%四次多项式拟合 y1=polyval(p1,x);y2=polyval(p2,x);%多项式求值plot(x,y,'c--',x,y1,'r:',x,y2,'y-.')p3=polyfit(x,y,2)%观察图像，类似抛物线，故用二次多项式拟合。

数值分析第九章计算实习题答案昆工(a)程序：clc;clear;a=1;b=2;%定义域h=0.05;%步长n=(b-a)/h;y0=1;%初值f= @(x,y) 1/x^2-y/x;%微分函数Xn=linspace(a,b,n+1);%将定义域分为n等份Yn=zeros(1,n);%结果矩阵Yn(1)=y0;%赋初值%以下根据改进欧拉公式求解for i=1:nxn=Xn(i);xnn=Xn(i+1);yn=Yn(i);yp=yn+h*f(xn,yn);yc=yn+h*f(xnn,yp);yn=(yp+yc)/2;Yn(i+1)=yn;endXn=Yn;%以下根据经典四阶R-K法公式求解for i=1:nxn=Xn(i);yn=Yn(i);k1=f(xn,yn);k2=f(xn+h/2,yn+h/2*k1);k3=f(xn+h/2,yn+h/2*k2);k4=f(xn+h,yn+h*k3);yn=yn+h/6*(k1+2*k2+2*k3+k4);Yn(i+1)=yn;enddisp(' 改进欧拉法四阶经典R-K法'); disp([Xn' Yn']) 结果如下：改进欧拉法四阶经典R-K法1 10.99887 0.998850.99577 0.99780.99114 0.996940.98532 0.996340.97857 0.996030.97111 0.996060.96311 0.996450.9547 0.997230.94598 0.998410.93705 10.92798 1.0020.91883 1.00440.90964 1.00730.90045 1.01060.89129 1.01430.88218 1.01840.87315 1.02290.86421 1.02780.85538 1.03310.84665 1.0388 （b）程序：clc;clear;a=0;b=1;%定义域H=[0.1 0.025 0.01];%步长y0=1/3;%初值f= @(x,y) -50*y+50*x^2+2*x;%微分函数xi=linspace(a,b,11);Y=1/3*exp(-50*xi)+xi.^2;%准确解Ym=zeros(1,11);for j=1:3h=H(j);n=(b-a)/h;Xn=linspace(a,b,n+1);%将定义域分为n等份Yn=zeros(1,n);%结果矩阵Yn(1)=y0;%赋初值for i=1:nxn=Xn(i);yn=Yn(i);k1=f(xn,yn);k2=f(xn+h/2,yn+h/2*k1);k3=f(xn+h/2,yn+h/2*k2);k4=f(xn+h,yn+h*k3);yn=yn+h/6*(k1+2*k2+2*k3+k4);Yn(i+1)=yn;endfor k=1:11m=0.1/h;Ym(k)=Yn(1+(k-1)*m);enddelta=Ym-Y;fprintf('步长为：%d \n', h);disp(' 四阶经典R-K法准确解误差'); disp([Ym' Y' delta']) end结果如下：步长为：1.000000e-01四阶经典R-K法准确解误差0.33333 0.33333 04.6055 0.012246 4.593263.062 0.040015 63.022864.05 0.09 863.9611844 0.16 118431.6235e+05 0.25 1.6235e+052.2256e+06 0.36 2.2256e+063.0509e+07 0.49 3.0509e+074.1823e+08 0.64 4.1823e+085.7333e+09 0.81 5.7333e+097.8594e+10 1 7.8594e+10步长为：2.500000e-02四阶经典R-K法准确解误差0.33333 0.33333 00.013015 0.012246 0.000768940.040063 0.040015 4.82e-050.090037 0.09 3.6857e-050.16004 0.16 3.6723e-050.25004 0.25 3.6722e-050.36004 0.36 3.6722e-050.49004 0.49 3.6722e-050.64004 0.64 3.6722e-050.81004 0.81 3.6722e-051 1 3.6722e-05步长为：1.000000e-02四阶经典R-K法准确解误差0.33333 0.33333 00.012256 0.012246 9.5673e-060.040016 0.040015 7.8252e-070.090001 0.09 6.6347e-070.16 0.16 6.6226e-070.25 0.25 6.6225e-070.36 0.36 6.6225e-070.49 0.49 6.6225e-070.64 0.64 6.6225e-070.81 0.81 6.6225e-071 1 6.6225e-07 由结果可知，步长越小，结果越精确。

数值分析计算实习题一学号：：院系：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的矩阵CA 的带元素ij a =C 中的元素1,i j s j c -++ 程序中A[5][501]即为压缩后的矩阵。

1.1.2 程序中的B[5][501]为过渡矩阵，在幂法迭代、反幂法迭代以及LU 分解中均用矩阵B 来计算，计算之间对B 进行适当的赋值。

插值法1.下列数据点的插值x 0 1 4 9 16 25 36 49 64y 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:nl=sym(y1(i));for k=1:i-1l=l*(x-x1(k))/(x1(i)-x1(k));endfor k=i+1:n l=l*(x-x1(k))/(x1(i)-x1(k));endLs=Ls+l; endLs=simplify(Ls) % 为所求插值多项式Ls(x).输出结果为Ls = -/*xA2+95549/72072*x-1/00*xA8-2168879/0*xA4+19/0*xA7+657859/*xA3+33983/ 0*xA5-13003/00*xA6(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=s plin e(x1,y1,x2);p=po Iyfit(x2,y3,3); % 得到三次样条拟合函数S=p(1)+p(2)*x+ p(3)*x^2+p(4)*xA3 % 得到S(x) 输出结果为：S =/6464-2399/88*x+/1984*xA2+2656867/624*xA3⑶ 在区间［0,64］上，分别对这两种插值和标准函数作图,Plot(x2,sqrt(x2),'b',x2,y2,'r',x2,y3,'y')蓝色曲线为y="X函数曲线，红色曲线为拉格朗日插值函数曲线，黄色曲线为三次样条插值曲线可以看到蓝色曲线与黄色曲线几乎重合，因此在区间[0,64] 上三次样条插值更精确。

弟二草插值法3.卜列数据点的插值可以得到平方根函数的近似，在区间064］上作图。

(1〉用这9个点做8次多项式插值Q x)。

(2)用三次样条(第一边界条件)程岸求S(X)。

从得到结果石在［0.64］ 1：・哪个插值更粘确：在区间［0,1］ I：•两种插值哪个更精确?(1) 8次多项式插值：(1)8次多项式插值:首先建立新的M-file:输入如卜代码(此为拉格朗口插值的功能函数)并保存function f=Language(x,y,x0)%求Li知数据点的拉格朗Fl插值多项式%己知数据点的x坐标向量：x%已知数据点的y坐标向量:y%插值的x坐标：x0%求得的拉格朗H插值多项式或在X0处的插值：fsyms t;ifi(lcngth(x)=length(y))n=length(x);elsedisp(*x和y的维数不相等!)；return;end %检错tbr(i=l:n)i=y(i);fbr(j=1:i-l)l=l*(t-x(j))/(x(i)-x(j))；end;for(j=i-M:n)end;for(j=i+l:n) l=l*(t-x(j))/(x(i)-x(j))； end;simplify(f);if(i==n) if|nargin=3)f=subs(C't\xO);else f=collcct(f);f=vpa(f,6);endendend再建立新的M-file：输入：clear;x=[0 1 49 16 25 36 49 64];y=[0：l：8]；%计算拉格朗口基丞数%计算拉格朗ri插值函数%化简%计算插值点的曲数值%将插值多项式展开%将插值多项式的系数化成6位精度的小数f=Uinguage(x,y) 运行得到f=1.32574*1-381410*t A2+.604294e-1 *t A3+.222972e-3 *t A5-.542921 e-5*t A6+.671268e・7T7・.328063e・9T8・.498071 e-2*t A4 这就是8次多项式插值L s(x)= 1.32574怜.381410*t A2+.604294e-1 *t A3+.222972e-3 *t A5-.542921 e-5*t A6+.671268e-7*t A7-.328063e-9*t A8-.498071 e-2*t A4. (2)三次样条插值：建立新的M-filc：输入：clear;x=[0 I 49 1625 36 4964];尸[0:8];t=[0:0.1:64];Y=t.A(0.5);O=Language(x,y)f= 1,32574*t-.381410*t.A2+.604294e-1 *t.A3+.222972e-3*t.A5-.542921 e・5*(. W+.671268e-7*t.A7-.328063e-9*t.A8-.498071 e-2 *t.A4;S=interp l(x,y,t.'spline,);plol(x,y,o；(・YY.lf.'b'」S'g:');grid;运行程序得到如下图：从结果屮很明显可以看出在[0.64].上.三次样条插值更精确，儿乎与原函数帀合。

《数值分析》计算实习作业《一》北航第一题 设有501501⨯的矩阵⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡=501500499321a bc b a b cc b a b ccb a bc c b a b c b a 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 矩阵的特征值)501,,2,1( =i i λ满足||min ||,501150121i i s λλλλλ≤≤=<<<试求1. 5011，λλ和s λ的值2. 的与数4015011λλκλμ-+=k 最接近的特征值)39,,2,1( =K κλi3. 的（谱范数）条件数2)A (cond 和行列式A det 要求1. 算法的设计方案（A 的所有零元素都不能存储）2. 全部源程序（详细注释）。

变量为double ，精度-1210=ε，输出为e 型12位有效数字3. 特征值s 5011，，λλλ和)39,,2,1( =K κλi 以及A cond det ,)A (2的值4. 讨论迭代初始向量的选取对计算结果的影响，并说明原因解答：1. 算法设计对于s λ满足||min ||5011i i s λλ≤≤=，所以s λ是按模最小的特征值，直接运用反幂法可求得。

对于5011，λλ，一个是最大的特征值，一个是最小的特征值，不能确定两者的绝对值是否相等，因此必须首先假设||||5011λλ≠，然后运用幂法，看能否求得一个特征值，如果可以求得一个，证明A 是收敛的，求得的结果是正确的，然后对A 进行带原点平移的幂法，偏移量是前面求得的特征值，可以求得另一个特征值，最后比较这两个特征值，较大的特征值是501λ，较小的特征值就是1λ。

如果在假设的前提下，无法运用幂法求得按模最大的特征值，即此时A 不收敛，则需要将A 进行带原点平移的幂法，平移量可以选取1，再重复上述步骤即可求得两个特征值。

插值法1.下列数据点的插值x 0 1 4 9 16 25 36 49 64y 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:nl=sym(y1(i));for k=1:i-1l=l*(x-x1(k))/(x1(i)-x1(k));endfor k=i+1:nl=l*(x-x1(k))/(x1(i)-x1(k));endLs=Ls+l;endLs=simplify(Ls) %为所求插值多项式Ls(x).输出结果为Ls =-24221063/63504000*x^2+95549/72072*x-1/3048192000*x^8-2168879/43545600 0*x^4+19/283046400*x^7+657859/10886400*x^3+33983/152409600*x^5-13003/2395 008000*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-2020406080100可以看到蓝色曲线与黄色曲线几乎重合，因此在区间[0,64]上三次样条插值更精确。

《数值分析》计算实习题目第一题：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 <stdio.h>#include <math.h>void init_a();//初始化Adouble get_an_element(int,int);//取A 中的元素函数double powermethod(double);//原点平移的幂法double inversepowermethod(double);//原点平移的反幂法int presolve(double);//三角LU 分解int solve(double [],double []);//解方程组int max(int,int);int min(int,int);double (*u)[502]=new double[502][502];//上三角U 数组double (*l)[502]=new double[502][502];//单位下三角L 数组double a[6][502];//矩阵Aint main(){int i,k;double lambdat1,lambdat2,lambda1,lambda501,lambdas,mu[40],det;double lambda[40];init_a();//初始化Alambdat1=powermethod(0);lambdat2=powermethod(lambdat1);lambda1=lambdat1<lambdat2?lambdat1:lambdat2;lambda501=lambdat1>lambdat2?lambdat1:lambdat2;presolve(0);lambdas=inversepowermethod(0);det=1;for(i=1;i<=501;i++)det=det*u[i][i];for (k=1;k<=39;k++){mu[k]=lambda1+k*(lambda501-lambda1)/40;presolve(mu[k]);lambda[k]=inversepowermethod(mu[k]);}printf("------------所有特征值如下------------\n");printf("λ=%1.11e λ=%1.11e\n",lambda1,lambda501);printf("λs=%1.11e\n",lambdas);printf("cond(A)=%1.11e\n",fabs(lambdat1/lambdas));printf("detA=%1.11e \n",det);for (k=1;k<=39;k++){printf("λi%d=%1.11e ",k,lambda[k]);if(k % 3==0) printf("\n");} delete []u;delete []l;//释放堆内存return 0;}void init_a()//初始化A{int i;for (i=3;i<=501;i++) a[1][i]=a[5][502-i]=-0.064;for (i=2;i<=501;i++) a[2][i]=a[4][502-i]=0.16;for (i=1;i<=501;i++) a[3][i]=(1.64-0.024*i)*sin(0.2*i)-0.64*exp(0.1/i); }double get_an_element(int i,int j)//从A中节省存储量的提取元素方法{if (fabs(i-j)<=2) return a[i-j+3][j];else return 0;}double powermethod(double offset)//幂法{int i,x1;double u[502],y[502];double beta=0,prebeta=-1000,yita=0;for (i=1;i<=501;i++)u[i]=1,y[i]=0;//设置初始向量u[]for (int k=1;k<=10000;k++){yita=0;for (i=1;i<=501;i++) yita=sqrt(yita*yita+u[i]*u[i]);for (i=1;i<=501;i++) y[i]=u[i]/yita;for (x1=1;x1<=501;x1++){u[x1]=0;for (int x2=1;x2<=501;x2++)u[x1]=u[x1]+((x1==x2)?(get_an_element(x1,x2)-offset):get_an_element(x1,x2))*y[x2] ;}prebeta=beta;beta=0;for (i=1;i<=501;i++) beta=beta+ y[i]*u[i];if (fabs((prebeta-beta)/beta)<=1e-12) {printf("offset=%f lambda=%f err=%e k=%d\n",offset,(beta+offset),fabs((prebeta-beta)/beta),k);break;};//输出中间过程,包括偏移量,误差,迭代次数}return (beta+offset);}double inversepowermethod(double offset)//反幂法{int i;double u[502],y[502];double beta=0,prebeta=0,yita=0;for (i=1;i<=501;i++)u[i]=1,y[i]=0; //设置初始向量u[]for (int k=1;k<=10000;k++){yita=0;for (i=1;i<=501;i++) yita=sqrt(yita*yita+u[i]*u[i]);for (i=1;i<=501;i++) y[i]=u[i]/yita;solve(u,y);prebeta=beta;beta=0;for (i=1;i<=501;i++) beta=beta+ y[i]*u[i];beta=1/beta;if (fabs((prebeta-beta)/beta)<=1e-12) {printf("offset=%f lambda=%f err=%e k=%d\n",offset,(beta+offset),fabs((prebeta-beta)/beta),k);break;};//输出中间过程,包括偏移量,误差,迭代次数}return (beta+offset);int presolve(double offset)//三角LU分解{int i,k,j,t;double sum;for (k=1;k<=501;k++)for (j=1;j<=501;j++){u[k][j]=l[k][j]=0;if (k==j) l[k][j]=1;} //初始化LU矩阵for (k=1;k<=501;k++){for (j=k;j<=min(k+2,501);j++){sum=0;for (t=max(1,max(k-2,j-2)) ; t<=(k-1) ; t++)sum=sum+l[k][t]*u[t][j];u[k][j]=((k==j)?(get_an_element(k,j)-offset):get_an_element(k,j))-sum;}if (k==501) continue;for (i=k+1;i<=min(k+2,501);i++){sum=0;for (t=max(1,max(i-2,k-2));t<=(k-1);t++)sum=sum+l[i][t]*u[t][k];l[i][k]=(((i==k)?(get_an_element(i,k)-offset):get_an_element(i,k))-sum)/u[k][k];}}return 0;}int solve(double x[],double b[])//解方程组{int i,t;double y[502];double sum;y[1]=b[1];for (i=2;i<=501;i++){sum=0;for (t=max(1,i-2);t<=i-1;t++)sum=sum+l[i][t]*y[t];y[i]=b[i]-sum;}x[501]=y[501]/u[501][501];for (i=500;i>=1;i--){sum=0;for (t=i+1;t<=min(i+2,501);t++)sum=sum+u[i][t]*x[t];x[i]=(y[i]-sum)/u[i][i];}return 0;}int max(int x,int y){return (x>y?x:y);}int min(int x,int y){return (x<y?x:y);}3.计算结果结果如下图所示:部分中间结果：给出了偏移量(offset)，误差(err)，迭代次数(k)4.讨论迭代初始向量的选取对计算结果的影响,并说明原因使用u[i]=1（i=1,2,...,501）作为初始向量进行迭代，可得出以上结果。

数值分析计算实习题二学号：姓名：院系：2015年11月28日一、算法采用拟上三角化的双步位移的QR分解法计算特征值，并利用列主元素Guass消去法求实数特征值对应的特征向量，定义全局变量n=10。

1 主程序在定义一维和二维数组时将数组的大小定为n+1（大小为11），这样在使用单个元素时的角标可以与矩阵或向量元素的角标相对应。

先输入矩阵A的各个元素；保留数组A，定义数组B的元素与A相同，即B=A，后面的矩阵变换用数组B来进行；对B进行拟上三角化；输出拟上三角化后的矩阵B；对拟上三角化后的B进行QR分解，输出矩阵Q、R、RQ；对拟上三角化的数组矩阵B进行QR分解；输出A的全部特征值；判断所有特征值，如果特征值为实数（即储存特征值虚部的数组中元素为0），利用列主元素高斯消去法求其特征向量；输出实特征值对应的特征向量。

2 拟上三角化程序对于r=1,2，……，n-2循环执行（1）求第r列的第r+2个元素到第n个元素的平方和，如果和为0，则跳出该次循环直接进行下一次r的循环。

如果和不为0，则继续执行。

（2）计算第r列从第r+1到第n个元素的平方和的平方根dr=当B[r+1][r]小于等于0时，取cr=dr，否则取cr=-drhr=cr*cr-cr*B[r+1][r]（3）令(0,...,0,[1][],[2][],...,[][])Tur B r r cr B r r B n r=+-+（4）计算()/()()/()()()/()()()()()()()T T Tpr B ur hr qr B ur hr tr pr ur hr wr qr tr ur B B wr ur ur pr ====-=--（5）执行下一次循环。

3 拟上三角化后的矩阵的QR 分解初始令Q=I 。

对于r=1,2，……，n-1循环执行（1）求第r 列的第r+1个元素到第n 个元素的平方和，如果和为0，则跳出该次循环直接进行下一次r 的循环。

第四章：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是积分下、上限%m是龙贝格积分表中行最大数目%wucha是两次相邻迭代值的绝对误差限n=1;h=b-a;wugu=1;x=a;k=0;RT=zeros(4,4);RT(1,1)=h*(feval(fun,a)+feval(fun,b))/2;while((wugu>wucha)&(k<m)|(k<4))k=k+1;h=h/2;s=0;for j=1:nx=a+h*(2*j-1);s=s+feval(fun,x);endRT(k+1,1)=RT(k,1)/2+h*s;n=2*n;for i=1:kRT(k+1,i+1)=((4^i)*RT(k+1,i)-RT(k,i))/(4^i-1);endwugu=abs(RT(k+1,k)-RT(k+1,k+1));endR=RT(k+1,k+1);输入：>>fun=inline('sqrt(x).*log(x)');>> [RT,R,wugu,h]=Romberg(fun,eps,1,1e-5,13)输出：RT =1 至5 列-0.000000268546145 0 0 0-0.245064670140209 -0.326752804004897 0 0-0.358104125949240 -0.395783944552250 -0.400386020588741 0 0-0.408090073087781 -0.424752055467295 -0.426683262861631 -0.427100679405645 0-0.429474601629505 -0.436602777810080 -0.437392825966266 -0.437562819031419 -0.437603847029951-0.438389494461832 -0.441361125405941 -0.441678348578999 -0.441746372747455 -0.4417627788404596 列-0.441766844267449R =-0.441766844267449wugu =4.065426989774412e-06h =0.031250000000000（3）自适应辛普森输入：>> f=inline('sqrt(x).*log(x)');>> q=quad(f,0,1,1e-4)输出：q =-0.4439755729517282.（1）复合辛普森建立m文件function q=combinesimpson2(F,x0,a,b,n)%复合Simpson多元求积公式%F—被积函数%x0—被积函数自变量%[a,b]积分区间%n—区间份数x=linspace(a,b,n+1);q=0;for k=1:nq=q+subs(F,x0,x(k))+4*subs(F,x0,(x(k)+x(k+1))/2)+subs(F,x0,x(k+1)); endq=q*(b-a)/n/6;输入：>> clear>> syms x y;>> F=exp(-x.*y);>> s=combinesimpson2(combinesimpson2(F,'x',0,1,4),'y',0,1,4)输出：s =exp(-1)/576 + exp(-1/2)/144 + exp(-1/4)/72 + exp(-3/4)/144 + exp(-1/8)/36 +exp(-3/8)/36 + exp(-5/8)/72 + exp(-7/8)/72 + (5*exp(-1/16))/144 + exp(-3/16)/24 + exp(-5/16)/36 + exp(-7/16)/36 + exp(-9/16)/144 + exp(-1/32)/36 + exp(-3/32)/18 + exp(-5/32)/36 + exp(-7/32)/36 + exp(-9/32)/36 + exp(-15/32)/36 + exp(-21/32)/36 + exp(-1/64)/36 + exp(-3/64)/18 + exp(-5/64)/18 + exp(-7/64)/18 + exp(-9/64)/36 + exp(-15/64)/18 + exp(-21/64)/18 + exp(-25/64)/36 + exp(-35/64)/18 + exp(-49/64)/36 + 47/576>> double(s)ans =0.796599967946203高斯求积公式function q=gaussquad(F,x0,a,b,n)%Gauss求积公式%F—被积函数%x0—被积函数自变量%[a,b]积分区间%n—节点个数syms t;F=subs(F,x0,(b-a)/2*t+(a+b)/2);[x,A]=gausspoints(n);q=(b-a)/2*sum(A.*subs(F,t,x));输入：>> clear>> syms x y;F=exp(-x.*y);>> s=gaussquad(gaussquad(F,x,0,1,4),y,0,1,4)输出：s =0.7966（2）复合辛普森输入：>> syms x y;>> f=exp(-x.*y);>> s=combinesimpson2(combinesimpson2(f,y,0,sqrt(1-x^2),4),x,0,1,4)输出：s =(3^(1/2)*(exp(-3^(1/2)/4) + 2*exp(-3^(1/2)/8) + 2*exp(-3^(1/2)/16) + 2*exp(-(3*3^(1/2))/16) + 4*exp(-3^(1/2)/32) + 4*exp(-(3*3^(1/2))/32) + 4*exp(-(5*3^(1/2))/32) + 4*exp(-(7*3^(1/2))/32) + 1))/576 + (7^(1/2)*(exp(-(3*7^(1/2))/16) + 2*exp(-(3*7^(1/2))/32) + 2*exp(-(3*7^(1/2))/64) + 2*exp(-(9*7^(1/2))/64) + 4*exp(-(3*7^(1/2))/128) + 4*exp(-(9*7^(1/2))/128) + 4*exp(-(15*7^(1/2))/128) + 4*exp(-(21*7^(1/2))/128) + 1))/1152 + (15^(1/2)*(exp(-15^(1/2)/16) + 2*exp(-15^(1/2)/32) + 2*exp(-15^(1/2)/64) + 2*exp(-(3*15^(1/2))/64) + 4*exp(-15^(1/2)/128) + 4*exp(-(3*15^(1/2))/128) + 4*exp(-(5*15^(1/2))/128) + 4*exp(-(7*15^(1/2))/128) + 1))/1152 + (15^(1/2)*(exp(-(7*15^(1/2))/64) + 2*exp(-(7*15^(1/2))/128) + 2*exp(-(7*15^(1/2))/256) + 2*exp(-(21*15^(1/2))/256) + 4*exp(-(7*15^(1/2))/512) + 4*exp(-(21*15^(1/2))/512) + 4*exp(-(35*15^(1/2))/512) + 4*exp(-(49*15^(1/2))/512) + 1))/1152 + (39^(1/2)*(exp(-(5*39^(1/2))/64) + 2*exp(-(5*39^(1/2))/128) + 2*exp(-(5*39^(1/2))/256) + 2*exp(-(15*39^(1/2))/256) + 4*exp(-(5*39^(1/2))/512) + 4*exp(-(15*39^(1/2))/512) + 4*exp(-(25*39^(1/2))/512) + 4*exp(-(35*39^(1/2))/512) + 1))/1152 + (55^(1/2)*(exp(-(3*55^(1/2))/64) + 2*exp(-(3*55^(1/2))/128) + 2*exp(-(3*55^(1/2))/256) + 2*exp(-(9*55^(1/2))/256) + 4*exp(-(3*55^(1/2))/512) + 4*exp(-(9*55^(1/2))/512) + 4*exp(-(15*55^(1/2))/512) + 4*exp(-(21*55^(1/2))/512) + 1))/1152 + (63^(1/2)*(exp(-63^(1/2)/64) + 2*exp(-63^(1/2)/128) + 2*exp(-63^(1/2)/256) + 2*exp(-(3*63^(1/2))/256) + 4*exp(-63^(1/2)/512) + 4*exp(-(3*63^(1/2))/512) + 4*exp(-(5*63^(1/2))/512) + 4*exp(-(7*63^(1/2))/512) + 1))/1152 + 1/24>> double(s)ans =0.670113633359095。

数值分析计算实习题第三章第二次作业：题一：x=-1:0.2:1;y=1./(1+25.*x.^2);f1=polyfit(x,y,3)f=poly2sym(f1)y1=polyval(f1,x)x2=linspace(-1,1,10)y2=interp1(x,y,x2)plot(x,y,'r*-',x,y1,'b-')hold onplot(x2,y2,'k')legend('数据点','3次拟合曲线','3次多项式插值')xlabel('X'),ylabel('Y')输出：f1 =0.0000 -0.5752 0.0000 0.4841f =(4591875547102675*x^3)/81129638414606681695789005144064 - (3305*x^2)/5746 + (1469057404776431*x)/20282409603651670423947251286016 + 4360609662300613/9007199254740992y1 =-0.0911 0.1160 0.2771 0.3921 0.4611 0.4841 0.4611 0.3921 0.2771 0.1160 -0.0911x2 =-1.0000 -0.7778 -0.5556 -0.3333 -0.1111 0.1111 0.3333 0.5556 0.7778 1.0000y2 =0.0385 0.0634 0.1222 0.3000 0.7222 0.7222 0.3000 0.1222 0.0634 0.0385题二：X=[0.0 0.1 0.2 0.3 0.5 0.8 1.0];Y=[1.0 0.41 0.50 0.61 0.91 2.02 2.46];p1=polyfit(X,Y,3)p2=polyfit(X,Y,4)Y1=polyval(p1,X)Y2=polyval(p2,X)plot(X,Y,'r*',X,Y1,'b-.',X,Y2,'g--')p3=polyfit(X,Y,2)Y3=polyval(p3,X)f1=poly2sym(p1)f2=poly2sym(p2)f3=poly2sym(p3)plot(X,Y,'r*',X,Y1,'b-.',X,Y2,'g--',X,Y3,'m--')legend('数据点','3次多项式拟合','4次多项式拟合','2次多项式拟合') xlabel('X轴'),ylabel('Y轴')输出：p1 =-6.6221 12.8147 -4.6591 0.9266p2 =2.8853 -12.3348 16.2747 -5.2987 0.9427Y1 =0.9266 0.5822 0.4544 0.5034 0.9730 2.0103 2.4602Y2 =0.9427 0.5635 0.4399 0.5082 1.0005 1.9860 2.4692p3 =3.1316 -1.2400 0.7356Y3 =0.7356 0.6429 0.6128 0.6454 0.8984 1.7477 2.6271f1 =- (7455778416425075*x^3)/1125899906842624 + (1803512222945435*x^2)/140737488355328 - (40981580032809*x)/8796093022208 + 8345953784399011/9007199254740992f2 =(1624271450198125*x^4)/562949953421312 - (3471944732519173*x^3)/281474976710656 + (4580931990070659*x^2)/281474976710656 - (1491459232922115*x)/281474976710656 + 1061409433081293/1125899906842624f3 =(18733*x^2)/5982 - (74179*x)/59820 + 73337/99700题三：建立三角插值函数的m文件function [A,B,Y1,Rm]=sanjiaobijin(X,Y,X1,m)%A B分别是m阶三角多项式Tm （x）的系数aj，bj（j=1,2,...,m）的系数矩阵，Y1是Tm（x）在X1处的值，X Y 数据点 ,Rm为均方误差n=length(X)-1;max1=fix((n-1)/2);if m>max1m=max1;endA=zeros(1,m+1);B=zeros(1,m+1);Ym=(Y(1)+Y(n+1))/2;Y(1)=Ym;Y(n+1)=Ym;A(1)=2*sum(Y)/n;for i=1:mB(i+1)=sin(i*X)*Y';A(i+1)=cos(i*X)*Y';endA=2*A/n;B=2*B/n;A(1)=A(1)/2;Y1=A(1);for k=1:mY1=Y1+A(k+1)*cos(k*X1)+B(k+1)*sin(k*X1);Tm=A(1)+A(k+1).*cos(k*X)+B(k+1).*sin(k*X);k=k+1;endY,Tm,Rm=(sum(Y-Tm).^2)/n输出：>> X=-pi:2*pi/33:pi;>> Y=X.^2.*cos(X);[A,B,Y1,Rm]=sanjiaobijin(X,Y,X1,16)输出：A =1 至 12 列-0.1397 4.4002 -2.8326 1.2355 -0.9128 0.7914 -0.7319 0.6982 -0.6773 0.6635 -0.6541 0.647413 至 17 列-0.6426 0.6393 -0.6370 0.6355 -0.6348B =1.0e-15 *1 至 12 列0 -0.0194 -0.0150 -0.0044 -0.0300 0.0105 0.0627 -0.0821 -0.0599 -0.0133 -0.0211 0.029713 至 17 列0.0178 0.0962 -0.1049 0.0328 -0.0122即可得16插值多项式的值X1=-pi:0.001:pi;[A,B,Y1,Rm]=sanjiaobijin(X,Y,X1,16)plot(X,Y,'r*',X1,Y1,'b-.')legend('数据点','16次三角插值多项式')xlabel('X轴'),ylabel('Y轴')。

1.下列数据点的插值x 0 1 4 9 16 25 36 49 64y 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:nl=sym(y1(i));for k=1:i-1l=l*(x-x1(k))/(x1(i)-x1(k));endfor k=i+1:nl=l*(x-x1(k))/(x1(i)-x1(k));endLs=Ls+l;endLs=simplify(Ls) %为所求插值多项式Ls(x).输出结果为Ls =-/*x^2+95549/72072*x-1/00*x^8-2168879/0*x^4+19/0*x^7+657859/*x^3+33983/ 0*x^5-13003/00*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 =/6464-2399/88*x+/1984*x^2+2656867/624*x^3(3)在区间[0,64]上，分别对这两种插值和标准函数作图，plot(x2,sqrt(x2),'b',x2,y2,'r',x2,y3,'y')蓝色曲线为y=函数曲线,红色曲线为拉格朗日插值函数曲线，黄色曲线为三次样条插值曲线可以看到蓝色曲线与黄色曲线几乎重合，因此在区间[0,64]上三次样条插值更精确。

数值分析计算实习题（二）数值分析计算实习题（二）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位。