数值分析作业
数值分析第一章作业
数值分析第一章作业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]区间上的分段线性插值多项式, 并写出误差估计式。
数值分析大作业
数值分析大作业数值分析大作业姓名：黄晨晨学号：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）求得结果相同，可验证结果正确。
数值分析作业
数值分析课后作业：习题一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.下列各近似值的绝对误差限都是1021⨯-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+bans =5 7 9>> b-aans =3 3 3>> a.*bans =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];>> cc = 1 2 3 1 2 3 1 2 3d = 4 5 6 4 5 6 4 5 6 >> c+dans =5 7 9 5 7 9 5 7 9>> d-cans = 3 3 33 3 33 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=1ax c1 x>0 其中 sgn = 0 x=0 -1 x<0 disp('Please input the coefficients of');disp('quadratic equation ax^2+bx+c=0, respectively') 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 Enddelta=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>0x 1=(-b-temp)/(2*a) End if b<0x 1=(-b+temp)/(2*a) End if b=0x 1=temp/(2*a) Endx 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)<epsdisp('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 Enddelta=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>0x 1=digit(digit(-b-temp,m)/digit(2*a,m),m) ; End if b<0x 1=digit(digit(-b+temp,m)/digit(2*a,m),m); End if b=0x 1=digit(temp/digit(2*a,m),m); Endx 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(11≤<--+∞=∑x nx nn n 和ln11...),12...53(2111253<<-++++++=-++x n x x x x x x n ，给出计算ln2的近似方法，编写相应的Matlab 程序，并比较算法运行情况. 解：方法一： x=1; s=0;for k=1:100s=s+(-1)^(k+1)*(x^k)/k; end sq=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所以方法二较方法一好。
数值分析作业
数值分析作业及答案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.用如下数值表构造不超过3次的插值多项式2. P55 11题.给出概率积分⎰-=xxdxey 022π的数据表用2次插值计算，试问：(1) 当x = 0.472时，积分值等于多少? (2) 当x 为何值时，积分值等于0.5？ 解：(1) 取x 0 = 0.47， x 1 = 0.48, x 2 = 0.4980.4955530040.04093346-80.1809899240.355496540.51166830.50274980.4937452=+=----⨯+----⨯+----⨯==----+----+----≈)48.049.0)(47.049.0()48.0472.0)(47.0472.0()49.048.0)(47.048.0()49.0472.0)(47.0472.0()49.0472.0)(48.047.0()49.0472.0)(48.0472.0()472.0())(())(())(())(())(())(()472.0(2120210221120121210y Lxx xx xxy x x x x x x y x x x x x x x x x x x x y(2)90.4769359350.05272367-80.4362204360.093439170.50274980.51166830.49374520.51166830.50274980.49374520.49 0.51166830.50274980.49374520.50274980.51166830.49374520.48 0.51166830.49374520.50274980.49374520.51166830.50274980.47=+=----⨯+----⨯+----⨯==----+----+----≈))(()5.0)(5.0())(()5.0)(5.0())(()5.0)(5.0()5.0())(())(())(())(())(())(()5.0(212210221120121210Lyyy y yy xy y yy y y xyyyy yy xy y y y y y x3. 证明方程e x +10x －2＝0在区间[0，1]内有一个根，如果使用二分法求该区间内的根，且误差不超过10－6，试问需要二分区间[0，1]多少次？4. 设x t ＝451.01为准确值，x a ＝451.023为x t 的近似值，试求出x a 有效数字的位数及相对误差 作业答案1.解：N 2(x ) = f (0)+f [0，1](x -0)+ f [0，1，2](x -0) (x -1) 1＋1×(x -0) ＋3×(x -0) (x -1)＝3x 2-2x +1 为求得P 3(x )，根据插值条件知，P 3(x )应具有下面的形式 P 3(x )＝N 2(x )＋k (x -0) (x -1) (x -2)，这样的P 3(x )自然满足：P 3(x i )= f (x i )由P 3’(1 )=3P 3’(1 )= N 2’(1 )+k (1-0) (1-2) =N 2’(1 )－k = 4－k=3∴ k =1∴ P 3(x )＝N 2(x )＋ (x -0) (x -1) (x -2)=x 3+1 3. 证明 令f (x )＝e x +10x －2，∵ f (0)=-1<0，f (1)=e+8> 0∴ f (x )= e x +10x －2 =0在[0，1]有根。
北理工数值分析大作业
数值分析上机作业第 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。
北航数值分析作业第一题
数值分析作业第一题一、 算法设计方案利用带状Dollittle 分解，将A[501][501]转存到数组C[5][501]，以节省存储空间1、计算λ1和λ501首先使用幂法求出矩阵的按模最大的特征值λ0：如果λ0>0，则其必为按模最大值，因此λ501=λ0，然后采用原点平移法，平移量为λ501，使用幂法迭代求出矩阵A -λ501I 的按模最大的特征值，其特征值按从小到大排列应为λ1-λ501、λ2-λ501、……、0。
因此A-λ501I 的按模最大的特征值应为λ1-λ501。
又因为λ501的值已求得，由此可直接求出λ1。
2、计算λSλS 为矩阵A 按模最小的特征值，可以通过反幂法直接求出。
3、计算λikλik 是对矩阵A 进行λik 平移后，再用反幂法求出按模最小的特征值λmin ，λik =λik +λmin 。
4、计算矩阵A 的条件数计算cond （A ）2和行列式det(A)矩阵A 的条件数为n12cond λλ）（ A ，其中λ1和λn 分别是矩阵A 的模最大和最小特征值，直接利用上面求得的结果直接计算。
矩阵A 的行列式可先对矩阵A 进行LU 分解后，det(A)等于U 所有对角线上元素的乘积。
二、源程序：#include<math.h>#include<stdio.h>#include<stdlib.h>#include<iostream.h>#define s 2#define r 2int Max(int v1,int v2);int Min(int v1,int v2);int maxt(int v1,int v2,int v3);void storage(double C[5][501],double b,double c);double mifa(double C[5][501]);void LU(double C[5][501]);double fmifa(double C[5][501]);int Max(int v1,int v2) //求两个数的最大值{ return((v1>v2)?v1:v2);}int Min(int v1,int v2) //求两个数最小值{ return ((v1<v2)?v1:v2);}int maxt(int v1,int v2,int v3) //求三个数最大值{ int t;if(v1>v2) t=v1;else t=v2;if(t<v3) t=v3;return(t);}/***将矩阵值转存在一个数组里，以节省存储空间***/void storage(double C[5][501],double b,double c){ int i=0,j=0;C[i][j]=0,C[i][j+1]=0;for(j=2;j<=500;j++)C[i][j]=c;i++;j=0;C[i][j]=0;for(j=1;j<=500;j++)C[i][j]=b;i++;for(j=0;j<=500;j++)C[i][j]=(1.64-0.024*(j+1))*sin(0.2*(j+1))-0.64*exp(0.1/(j+1));i++;for(j=0;j<=499;j++)C[i][j]=b;C[i][j]=0;i++;for(j=0;j<=498;j++)C[i][j]=c;C[i][j]=0,C[i][j+1]=0;}//用于求解最大的特征值，幂法double mifa(double C[5][501]){ int m=0,i,j;double b2,b1=0,sum;double u[501],y[501];for (i=0;i<501;i++){ u[i] = 1.0;}do{ sum=0;if(m!=0)b1=b2;m++;for(i=0;i<=500;i++)sum+=u[i]*u[i];for(i=0;i<=500;i++)y[i]=u[i]/sqrt(sum);for(i=0;i<=500;i++){ u[i]=0;for(j=Max(i-r,0);j<=Min(i+s,500);j++)u[i]=u[i]+C[i-j+s][j]*y[j];}b2=0;for(i=0;i<=500;i++)b2=b2+y[i]*u[i];}while(fabs(b2-b1)/fabs(b2)>=1.0e-12);return b2;}/*****行列式LU分解*****/void LU(double C[5][501]){ double sum;int k,i,j;for(k=1;k<=501;k++){ for(j=k;j<=Min(k+s,501);j++){ sum=0;for(i=maxt(1,k-r,j-s);i<=k-1;i++)sum+=C[k-i+s][i-1]*C[i-j+s][j-1];C[k-j+s][j-1]-=sum;}for(j=k+1;j<=Min(k+r,501);j++){ sum=0;for(i=maxt(1,j-r,k-s);i<=k-1;i++)sum+=C[j-i+s][i-1]*C[i-k+s][k-1];C[j-k+s][k-1]=(C[j-k+s][k-1]-sum)/C[s][k-1];}}}/***带状DOOLITE分解，并且求解出方程组的解***/void solve(double C[5][501],double x[501],double b[501]){ int i,j,k,t;double B[5][501],c[501];for(i=0;i<=4;i++){ for(j=0;j<=500;j++)B[i][j]=C[i][j];}for(i=0;i<=500;i++)c[i]=b[i];for(k=0;k<=500;k++){ for(j=k;j<=Min(k+s,500);j++){ for(t=Max(0,Max(k-r,j-s));t<=k-1;t++)B[k-j+s][j]=B[k-j+s][j]-B[k-t+s][t]*B[t-j+s][j];}for(i=k+1;i<=Min(k+r,500);i++){ for(t=Max(0,Max(i-r,k-s));t<=k-1;t++)B[i-k+s][k]=B[i-k+s][k]-B[i-t+s][t]*B[t-k+s][k];B[i-k+s][k]=B[i-k+s][k]/B[s][k];}}for(i=1;i<=500;i++)for(t=Max(0,i-r);t<=i-1;t++)c[i]=c[i]-B[i-t+s][t]*c[t];x[500]=c[500]/B[s][500];for(i=499;i>=0;i--){ x[i]=c[i];for(t=i+1;t<=Min(i+s,500);t++)x[i]=x[i]-B[i-t+s][t]*x[t];x[i]=x[i]/B[s][i];}}//用于求解模最大的特征值，反幂法double fmifa(double C[5][501]){ int m=0,i;double b2,b1=0,sum=0,u[501],y[501];for (i=0;i<=500;i++){ [i] = 1.0;}do{ if(m!=0)b1=b2;m++;sum=0;for(i=0;i<=500;i++)sum+=u[i]*u[i];for(i=0;i<=500;i++)y[i]=u[i]/sqrt(sum);solve(C,u,y);b2=0;for(i=0;i<=500;i++)b2+=y[i]*u[i];}while(fabs(b2-b1)/fabs(b2)>=1.0e-12);return 1/b2;}/***主程序***/void main(){ double b=0.16,c=-0.064,det=1.0;int i;double C[5][501],cond;storage(C,b,c); //进行C的赋值cout.precision(12); //定义输出精度double k1=mifa(C); //利用幂法计算矩阵的最大特征值和最小特征值if(k1<0)printf("λ1=%.12e\n",k1);else if(k1>=0)printf("λ501=%.12e\n",k1);for(i=0;i<501;i++)C[2][i]=C[2][i]-k1;double k2=mifa(C)+k1;if(k2<0)printf("λ1=%.12e\n",k2);else if(k2>=0)printf("λ501=%.12e\n",k2);storage(C,b,c);double k3=fmifa(C); //利用反幂法计算矩阵A的按模最小特征值printf("λs=%.12e\n",k3);storage(C,b,c); //计算最接近特征值double u[39]={0};for(i=0;i<39;i++){ u[i]=k1+(i+1)*(k2-k1)/40;C[2][i]=C[2][i]-u[i];u[i]=fmifa(C)+u[i];printf("与数u%d 最接近的特征值λ%d: %.12e\n",i+1,i+1,u[i]);}if(k1>0) //计算矩阵A的条件数，取2范数cond=fabs(k1/k3);else if(k1<0)cond=fabs(k2/k3);storage(C,b,c);LU(C); //利用LU分解计算矩阵A的行列式for(i=0;i<501;i++)det*=C[2][i];printf("\ncond(A)=%.12e\n",cond);printf("\ndet(A)=%.12e\n",det);}三、计算结果：四、结果分析迭代初始向量的选择对果有一定的影响，选择不同的初始向量可能会得到不同阶的特征值。
数值分析作业答案——第五版
数值分析第一次作业及参考答案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、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 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)