数值分析作业
数值分析第一章作业
数值分析第一章作业1.数值计算方法设计的基本手段是( ).(A) 近似 (B) 插值 (C) 拟合 (D) 迭代2.为了在有限时间内得到结果,用有限过程取代无限过程所产生的近似解与精确解之间的误差称为( ).(A) 舍入误差 (B) 截断误差 (C) 测量误差 (D) 绝对误差3.由于计算机的字长有限,原始数据在机器内的表示以及进行算术运算所产生的误差统称为( ).(A) 舍入误差 (B) 截断误差 (C) 相对误差 (D) 绝对误差4.数值计算方法研究的核心问题可以概括为( )对计算结果的影响.(A) 算法的稳定性 (B) 算法的收敛性 (C) 算法的复杂性 (D) 近似5.当N 充分大时,利用下列各式计算121N N dx I x+=+⎰,等式( )得到的结果最好. (A) arctan(1)arctan()I N N =+- (B) 2arctan(1)I N N =++ (C) 21arctan()1I N N =++ (D) 211I N =+6.计算61), 1.4≈,利用下列哪个公式得到的结果最好?为什么?(B) 3(3- (D) 99-7.计算球体的体积,已知半径的相对误差限不超过3310-⨯,则计算所得体积的相对误差限如何估计?8.设0x >,近似值*x 的相对误差限为δ,试估计*ln x 的误差限.9.计算圆柱体的体积,已知底面半径r 及圆柱高h 的相对误差限不超过δ,则计算所得体积的相对误差限如何估计?.10.用秦九韶算法求32()431f x x x x =-+-在2x =处的值.11.已知近似值 1.0000x *=的误差限4()110x ε*-=⨯,21()16f x x =,求(())f x ε*,并说明x *及()f x *的各有几位有效数字.12.设a 为非零常数,已知0y 的近似值0y *,由递推式1n n y ay -=计算序列{}n y 的近似值,分析该算法的稳定性.。
数值分析大作业一
数值分析大作业一一、算法设计方案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代表端点的一阶导。
清华大学高等数值分析 第一次实验作业
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
数值分析第二次作业
数值分析大作业（二）院系：学号：姓名：一、算法设计方案：1、矩阵A的存储和特征值的数据结构设计A是一个10×10的方阵，且由给定的元素计算公式知A是非对称的矩阵，因此采取二维数组来存储数组A。
考虑到特征值为复数的情况，采用结构体来存储特征值的实部和虚部，当特征值为实数时，虚部赋值为零。
2、求解过程的分析首先利用HouseHolder矩阵将A拟上三角化，得到A(n-1)。
然后利用带双步位移的QR分解法求解矩阵A的特征值。
利用选列主元素的高斯消去法求解A的实特征值对应的特征向量，将等式两边移项后，相当于求解齐次线性方程组。
在计算前首先要求齐次方程组的系数矩阵的秩，从而确定自由量,赋予适当的初值，从而可以求得特征值对应的一个特征值。
利用QR分解的基本方法可以得到A(n-1)进行QR分解后的Q、R和RQ。
二、源程序#include <stdio.h>#include <math.h>#define N 10#define epsilon 1E-12#define L 10000typedef struct Comp{double real;double image;}Comp;void getRealMatrix(int n,double A[][n],double Q[][n],double regulate); int isContainThisValue(int m,int mainIndex[],int index);int maxElementIndex(int m,int n,int colNum,double array[][n],int mainIndex[]);void mainElementElimination(int n,double equationArray[][n],double solveArray[],double);void qrMethod(int,double a[][]);int sgn(double realNum);void outputMatrix(int m,int n,double A[][n]);void calcMatrix(int m,int n,double A[][n],double M[][n],double s,double t);void productMatrix(int m,int n,double A[][n],double M[][n]);void hessenbergMatrix(int n,double matrix_a[][]);void solvePowerEquation(double b,double c,Comp twoDArr[2]);void qrIterateFormula(int m,int n,double M[][n],double A[][n]);void initialMatrix(int n, double A[][]);double crossProduct(double *p_1,double *p_2,int startIndex,int endIndex,int increment);void transpositionMatirx(int m,int n,double matrix[][]);void iterateFormula(int m,int n,double matrix[][],double w[],double u[],double p[]);void qrWithTwoStepMove(int n,double matrix[][],double matrix_Q[][],Comp lamdaValue[],double e,double iterCount);int main(){double Matrix_Q[N][N],Matrix_A[N][N];Comp lamdas[N];double vector[N]={0};int i,j;//初始化矩阵initialMatrix(N,Matrix_A);//拟上三角化化hessenbergMatrix(N,Matrix_A);//输出拟上三角化后的矩阵printf("A拟上三角化的矩阵：\n");outputMatrix(N,N,Matrix_A);//双步位移的QRqrWithTwoStepMove(N,Matrix_A,Matrix_Q,lamdas,epsilon,L);//求解Q、R、RQinitialMatrix(N,Matrix_A);hessenbergMatrix(N,Matrix_A);qrMethod(N,Matrix_A);for(i=0;i<=N-1;i++){printf("%.12E,%.12E\n",lamdas[i].real,lamdas[i].image);}//求解实特征值对应的特征向量initialMatrix(N,Matrix_A);for(i=0;i<=N-1;i++){if(lamdas[i].image!=0)continue;else{getRealMatrix(N,Matrix_A,Matrix_Q,lamdas[i].real);mainElementElimination(N,Matrix_Q,vector,1);printf("%.12E的特征向量为\n",lamdas[i].real);for(j=0;j<=N-1;j++)printf("%.12E ",vector[j]);printf("\n");}}return 0;}//输出矩阵void outputMatrix(int m,int n,double A[][n]){int i,j;for(i=0;i<=m-1;i++){for(j=0;j<=m-1;j++){printf("%0.12E ",A[i][j]);}printf("\n");}printf("\n");}//初始化矩阵的元素void initialMatrix(int n, double A[][n]){int i,j;for(i=0;i<=n-1;i++)for(j=0;j<=n-1;j++){if(i==j)A[i][j] = 1.5*cos((i+1)+1.2*(j+1));elseA[i][j] = sin(0.5*(i+1)+0.2*(j+1));}}//对矩阵进行拟上三角化void hessenbergMatrix(int n,double matrix_a[][n]){int r,i;double u[n],p[n],w[n],c,d,h;for(r=0;r<=n-3;r++){if(crossProduct(&matrix_a[r+2][r],&matrix_a[r+2][r],r+2,n-1,n) == 0)continue;d = sqrt(crossProduct(&matrix_a[r+1][r],&matrix_a[r+1][r],r+1,n-1,n)); c = matrix_a[r+1][r] == 0?d:(-sgn(matrix_a[r+1][r])*d);h = c*c - c*matrix_a[r+1][r];for(i=0;i<=n-1;i++){u[i] = (i<(r+1))?0:(i==(r+1))?(matrix_a[r+1][r]-c):matrix_a[i][r];}transpositionMatirx(n,n,matrix_a);for(i=0;i<=n-1;i++){p[i] = crossProduct(&matrix_a[0][0]+n*i,u,0,n-1,1)/h;}transpositionMatirx(n,n,matrix_a);for(i=0;i<=n-1;i++){w[i] = (crossProduct(&matrix_a[0][0]+n*i,u,0,n-1,1)- crossProduct(p,u,0,n-1,1)*u[i])/h;}iterateFormula(n,n,matrix_a,w,u,p);}}//双步位移的QR算法实现void qrWithTwoStepMove(int n,double matrix[][n],double matrix_Q[][n],Comp lamdaValue[],double e,double iterCount){int k,m,i=0,index=0;double s,t;Comp solver[2];k =1;m = n-1;while(1){if(fabs(matrix[m][m-1])<=e){lamdaValue[index].real = matrix[m][m];lamdaValue[index].image =0;index++;m=m-1;if(m<2){if(m==0){lamdaValue[index].real= matrix[m][m];lamdaValue[index].image = 0;index++;}else if(m==1){s = matrix[m-1][m-1]+ matrix[m][m];t = matrix[m-1][m-1]*matrix[m][m] - matrix[m][m-1]*matrix[m-1][m];solvePowerEquation(s,t,solver);for(i=0;i<=1;i++,index++)lamdaValue[index] = solver[i];}return;}}else{if(fabs(matrix[m-1][m-2])<=e){s = matrix[m-1][m-1]+ matrix[m][m];t = matrix[m-1][m-1]*matrix[m][m] - matrix[m][m-1]*matrix[m-1][m];solvePowerEquation(s,t,solver);for(i=0;i<=1;i++,index++)lamdaValue[index] = solver[i];m=m-2;if(m<2){if(m==0){lamdaValue[index].real= matrix[m][m];lamdaValue[index].image = 0;index++;}else if(m==1){s = matrix[m-1][m-1]+ matrix[m][m];t = matrix[m-1][m-1]*matrix[m][m] - matrix[m][m-1]*matrix[m-1][m];solvePowerEquation(s,t,solver);for(i=0;i<=1;i++,index++)lamdaValue[index] = solver[i];}return;}}else{if(k == iterCount)return;else{s = matrix[m-1][m-1]+ matrix[m][m];t = matrix[m-1][m-1]*matrix[m][m] - matrix[m][m-1]*matrix[m-1][m];calcMatrix(m,n,matrix,matrix_Q,s,t);qrIterateFormula(m+1,n,matrix_Q,matrix);k =k + 1;}}}}}//矩阵拟上三角化的迭代式子void iterateFormula(int m,int n,double matrix[][n],double w[],double u[],double p[]){int i,j;for(i=0;i<=m-1;i++)for(j=0;j<=m-1;j++){matrix[i][j] = matrix[i][j]-w[i]*u[j] - u[i]*p[j];}}//计算向量的内积， startIndex是起始下标，endIndex是结束下标，increment是元素增量double crossProduct(double *p_1,double *p_2,int startIndex,int endIndex,int increment){double sum = 0;int i;if(p_1 == p_2){for(i=startIndex;i<=endIndex;i++){sum += (*p_1)*(*p_1);p_1 = p_1 + increment;}}else{for(i=startIndex;i<=endIndex;i++){sum += (*p_1)*(*p_2);p_1 = p_1 + increment;p_2 = p_2 + increment;}}return sum;}//矩阵转置void transpositionMatirx(int m,int n,double matrix[][n]){int i,j;double changeVar;for(i=0;i<=m-1;i++){for(j=i;j<=m-1;j++){if(i != j){changeVar = matrix[i][j];matrix[i][j] = matrix[j][i];matrix[j][i] = changeVar;}}}}//符号函数int sgn(double realNum){return realNum== 0? 0 : (realNum>0?1:-1);}//代入的矩阵void calcMatrix(int m,int n,double A[][n],double M[][n],double s,double t){int i,j;productMatrix(m,n,A,M);for(i=0;i<=m;i++)for(j=0;j<=m;j++){M[i][j] = M[i][j] - s*A[i][j] + ((i==j)? t:0);}}//计算A平方void productMatrix(int m,int n,double A[][n],double M[][n]){int i,j,r;double sum;for(i=0;i<=m;i++)for(j=0;j<=m;j++){sum = 0;for(r=0;r<=m;r++)sum += A[i][r]*A[r][j];M[i][j]= sum;}}//计算lamdavoid solvePowerEquation(double b,double c,Comp twoDArr[2]){double delta = b*b-4*c;if(delta<0){twoDArr[0].real = b/2;twoDArr[0].image = sqrt(fabs(delta))/2;twoDArr[1].real = b/2;twoDArr[1].image = -sqrt(fabs(delta))/2;}else{twoDArr[0].real = (b+sqrt(delta))/2;twoDArr[0].image = 0;twoDArr[1].real = (b-sqrt(delta))/2;twoDArr[1].image = 0;}}//双步移位QR中的迭代公式void qrIterateFormula(int m,int n,double M[][n],double A[][n]) {int r,i;double c,d,h,u[m],v[m],p[m],w[m];for(r=0;r<=m-2;r++){if(crossProduct(&(M[r+1][r]),&(M[r+1][r]),r+1,m-1,n) == 0) {continue;}else{d = sqrt(crossProduct(&M[r][r],&M[r][r],r,m-1,n));c = (M[r][r] == 0)?d:(-sgn(M[r][r])*d);h = c*c - c*M[r][r];for(i=0;i<=m-1;i++){u[i] = (i<r)?0:(i==r)?(M[r][r]-c):M[i][r];}transpositionMatirx(m,n,M);for(i=0;i<=m-1;i++){v[i] = crossProduct(&M[0][0]+n*i,u,0,m-1,1)/h;}transpositionMatirx(m,n,M);for(i=0;i<=m-1;i++)p[i]=0;iterateFormula(m,n,M,u,v,p);//第二部分计算transpositionMatirx(m,n,A);for(i=0;i<=m-1;i++){p[i] = crossProduct(&A[0][0]+n*i,u,0,m-1,1)/h;}transpositionMatirx(m,n,A);for(i=0;i<=m-1;i++){w[i] = (crossProduct(&A[0][0]+n*i,u,0,m-1,1) - crossProduct(p,u,0,m-1,1)*u[i])/h;}iterateFormula(m,n,A,w,u,p);}}}//基本QR方法void qrMethod(int n,double A[][n]){int r,i,j;double c,d,h,sum;double u[n],w[n],p[n];double Q[n][n];double R[n][n];//初始化为单位阵for(i=0;i<=n-1;i++)for(j=0;j<=n-1;j++){if(i==j){Q[i][j]=1;R[i][j] = 0;}else{Q[i][j]=0;R[i][j] = 0;}}for(r=0;r<=n-1;r++){//判断第r列元素从 r+1行开始是否全部为0if(crossProduct(&A[r+1][r],&A[r+1][r],r+1,n-1,n)==0)continue;else{d = sqrt(crossProduct(&A[r][r],&A[r][r],r,n-1,n));c = (A[r][r] == 0)?d:(-sgn(A[r][r])*d);h = c*c - c*A[r][r];//初始化向量ufor(i=0;i<=n-1;i++){u[i] = (i<r)?0:(i==r)?(A[r][r]-c):A[i][r];}for(i=0;i<=n-1;i++){w[i] = crossProduct(&Q[0][0]+n*i,u,0,n-1,1)/h;}//辅助作用for(i=0;i<=n-1;i++)p[i]=0;//迭代出QiterateFormula(n,n,Q,w,u,p);transpositionMatirx(n,n,A);for(i=0;i<=n-1;i++){p[i] = crossProduct(&A[0][0]+n*i,u,0,n-1,1)/h;}//辅助作用for(i=0;i<=n-1;i++)w[i]=0;transpositionMatirx(n,n,A);iterateFormula(n,n,A,w,u,p);}}printf("正交阵Q: \n");outputMatrix(n,n,Q);printf("上三角阵R: \n");outputMatrix(n,n,A);for(i=0;i<=n-1;i++){for(j=0;j<=n-1;j++){sum=0;for(r=i;r<=n-1;r++)sum += A[i][r]*Q[r][j];R[i][j] = sum;}}printf("输出RQ：\n");outputMatrix(N,N,R);}//得到齐次线性方程组的系数矩阵void getRealMatrix(int n,double A[][n],double Q[][n],double regulate) {int i,j;for(i=0;i<=n-1;i++)for(j=0;j<=n-1;j++){if(i==j)Q[i][j] = A[i][j] - regulate;elseQ[i][j] = A[i][j];}}//列主元素消元法void mainElementElimination(int n,double equationArray[][n],double solveArray[],double referValue){int mainIndex[n];int temp;int k=0,i=0,j=0;double sum = 0;double ratio,afterDivide;for(i=0;i<=n-1;i++)mainIndex[i] = -1;for(k=0;k<=n-1;k++){//寻找最大元素所在的行号mainIndex[k]=maxElementIndex(n,n,k,equationArray,mainIndex);temp = mainIndex[k];for(i=0;i<=n-1;i++){if(isContainThisValue(n,mainIndex,i))continue;else{ratio = equationArray[i][k]/equationArray[mainIndex[k]][k];for(j=k;j<=n-1;j++){equationArray[i][j] = equationArray[i][j] - equationArray[mainIndex[k]][j]*ratio;afterDivide = equationArray[i][j];}}}}//代入一个设定值solveArray[n-1] = referValue;for(i=n-2;i>=0;i--){sum = 0;for(j=i+1;j<=n-1;j++){sum += equationArray[mainIndex[i]][j]*solveArray[j];}solveArray[i] = - sum/equationArray[mainIndex[i]][i];}}//返回一列中最大元素的下标int maxElementIndex(int m,int n,int colNum,double array[][n],int mainIndex[]){int maxIndex=0;int i=0,j=0;double max = array[0][colNum];for(i=0;i<=m-1;i++){if(!isContainThisValue(m,mainIndex,i)){max = array[i][colNum];maxIndex = i;break;}}for(i++;i<=m-1;i++){if(isContainThisValue(m,mainIndex,i))continue;if(max<array[i][colNum]){max = array[i][colNum];maxIndex = i;}}return maxIndex;}//判断index是否在mainIndex数组中int isContainThisValue(int m,int mainIndex[],int index) {int i = 0;for(i=0;i<=m-1;i++){if(index == mainIndex[i])return 1;}return 0;}三、计算结果。
数值分析作业题(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]区间上的分段线性插值多项式, 并写出误差估计式。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
return q;
}
double Rombf(double x)
{
double y;
y=1/x;
return y;
}
main()
{
double a=1;
double b=3;
double eps=1e-5;
double t=Romb(a,b,eps);
printf("The result is:%.3f\n",t);
k=13 p1=1.0000 err=2.6008e-005 y= 2.0293e-009
k=14 p1=1.0000 err=1.3004e-005 y=5.0732e-010
k=15 p1 =1.0000 err=6.5020e-006 y=1.2683e-010
k=16 p1 =1.0000 err=3.2510e-006 y=3.1708e-011
k=17 p1 =1.0000 err=1.6255e-006 y =7.9272e-012
k=18 p1 =1.0000 err =8.1279e-007 y= 1.9820e-012
ans = 1.0000
结果说明：经过18次迭代得到精确解为1，误差为8.1279e-007。
第三章
1.题目：
用最小二乘法求一个形如y=a+bx^2的经验公式,使他与下列数据相拟合,并求均方误差
3.实验结果
p0 = 1.2000
k =1 p1=1.1030 err=0.0970 y=0.0329
k= 2 p1=1.0524 err=0.0507 y=0.0084
k =3 p1=1.0264 err=0.0260 y=0.0021
k =4 p1=1.0133 err=0.0131 y=5.2963e-004
2.程序如下：
>> x=[19 25 31 38 44];
>> y=[19.0 32.3 49.0 73.3 97.8];
>> [p,v]=ZXRC_Poly(x.^2,y,1)
p =
0.05003512421916010.972578656906791
v =
0.122569200640555
3. 运行结果：
>> b=[0.4127;1.7321;-0.8621];
>> x=GaJo_inv(A)*b
x =4.58668603133057
-0.63152317337598
2.73520013957385
>>A*x-b
ans =
1.0e-015 *
0.05551115123126
0.22204460492503
k =5 p1=1.0066 err=0.0066 y=1.3270e-004
k =6 p1=1.0033 err=0.0033 y=3.3211e-005
k =7 p1=1.0017 err=0.0017 y=8.3074e-006
k =8 p1=1.0008 err=8.3157e-004 y = 2.0774e-006
第二章
1.题目：运用MATLAB编程实现牛顿迭代
2. 实验操作
1、打开MATLAB程序软件。
2、在MATLAB中编辑如下的M程序。
function [p1,err,k,y]=newton(f,df,p0,delta,max)
%f是要求根的方程(f(x)=0);
%df是f(x)的导数；
%p0是所给初值,位于x*附近；
1.57142857142857
2.14285714285714
2.71428571428571
2.28571428571429
1.85714285714286
1.42857142857143
第六章
1.题目：
解下列方程组
2.解答过程：
>>A=[0.6428 0.3475 -0.8468;0.3475 1.8423 0.4759;-0.8468 0.4759 1.2147];
{
x=a+(i+0.5)*h;
p=p+Rombf(x);
}
p=(y[0]+h*p)/2.0;
s=1.0;
for(k=1;k<=m;k++)
{
s=4.0*s;
q=(s*p-y[k-1])/(s-1.0);
y[k-1]=p;p=q;
}
ep=fabs(q-y[m-1]);
m=m+1;y[m-1]=q;n=n+n;h=h/2.0;
{
double Rombf();
int m,n,i,k;
double y[10],h,ep,p,x,s,q;
h=b-a;
y[0]=h*(Rombf(a)+Rombf(b))/2.0;
m=1;n=1;ep=eps+1.0;
while((ep>=eps)&&(m<=9))
{
p=0.0;
for(i=0;i<=n-1;i++)
%delta是给定允许误差；
%max是迭代的最大次数；
%p1是newton法求得的方程的近似解；
%err是p0的误差估计；
%k是迭代次数；
p0
for k=1:max
p1=p0-feval('f',p0)/feval('df',p0);
err=abs(p1-p0);
p0=p1;
k
p1
err
y=feval('f',p1)
if (err<delta)|(y==0)|(k==max)
break;
end
end
4、M文件2
function y=f(x)
y=x^3-3*x+2;
5、M文件3
function y=df(x)
y=3*x^2-3;
6、在程序窗口中，调用上面的M文件，对具体问题求解
>>newton('f','df',1.2,10^(-6),20)
k =9 p1=1.0004 err=4.1596e-004 y =5.1943e-007
k=10 p1=1.0002 err=2.0802e-004 y= 1.2987e-007
k=11 p1=1.0001 err=1.0402e-004 y =3.2468e-008
k=12 p1=1.0001 err=5.2014e-005 y=8.1170e-009
0.55511151231258
}
3．运行结果：
The result is:1.099
第五章
1.题目：用追赶法解三对角矩阵方程AX=b,其中:
2.程序如下：
>> a=-1*ones(1,5);
>> c=a;
>> b=2*ones(1,6);
>> f=[1 0 1 0 0 1];
>> [x]=ZhuiGan(a,b,c,f)’
x =
a=0.972578656906791 b=0.05003512421章
1.题目：用龙贝格法计算积分： ，要求
2.程序如下：
#include<stdio.h>
#include<math.h>
double Romb(double a,double b,double eps)