数值分析第五版计算实习题
数值分析计算实习题二
《数值分析》计算实习题二算法设计方案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，该函数可对系数矩阵相同// 而右端向量不同的多个方程组同时求解。
最新数值分析(第五版)计算实习题第五章作业
数值分析第五章第一题：LU分解法：建立m文件function h1=zhijieLU(A,b)%h1各阶主子式的行列式值[n n]=size(A);RA=rank(A);if RA~=ndisp('请注意：因为A的n阶行列式h1等于零，所以A不能进行LU分解。
A的秩RA如下：')RA,h1=det(A);returnendif RA==nfor p=1:nh(p)=det(A(1:p,1:p));endh1=h(1:n);for i=1:nif h(1,i)==0disp('请注意：因为A的r阶主子式等于零，所以A不能进行LU分解。
A的秩RA和各阶顺序主子式h1依次如下：')h1;RAreturnendendif h(1,i)~=0disp('请注意：因为A的r阶主子式都不等于零，所以A能进行LU分解。
A的秩RA和各阶顺序主子式h1依次如下：')for j=1:nU(1,j)=A(1,j);endfor k=2:nfor i=2:nfor j=2:nL(1,1)=1;L(i,i)=1;if i>jL(1,1)=1;L(2,1)=A(2,1)/U(1,1);L(i,1)=A(i,1)/U(1,1);L(i,k)=(A(i,k)-L(i,1:k-1)*U(1:k-1,k))/U(k,k);elseU(k,j)=A(k,j)-L(k,1:k-1)*U(1:k-1,j);endendendendh1;RA,U,L,X=inv(U)*inv(L)*bendend输入：>> A=[10 -7 0 1;-3 2.099999 6 2;5 -1 5 -1;2 1 0 2];>> b=[8;5.900001;5;1];>> h1=zhijieLU(A,b)输出：请注意：因为A的r阶主子式都不等于零，所以A能进行LU分解。
A的秩RA和各阶顺序主子式h1依次如下：RA =4U =10.0000 -7.0000 0 1.00000 2.1000 6.0000 2.30000 0 -2.1429 -4.23810 -0.0000 0 12.7333L =1.0000 0 0 0-0.3000 1.0000 0 00.5000 1.1905 1.0000 -0.00000.2000 1.1429 3.2000 1.0000X =-0.2749-1.32981.29691.4398h1 =10.0000 -0.0000 -150.0001 -762.0001列主元高斯消去法：建立m文件function [RA,RB,n,X]=liezhu(A,b)B=[A b];n=length(b);RA=rank(A);RB=rank(B);zhicha=RB-RA;if zhicha>0disp('请注意：因为RA~=RB，所以方程组无解')returnwarning off MATLAB:return_outside_of_loopendif RA==RBif RA==ndisp('请注意：因为RA=RB，所以方程组有唯一解')X=zeros(n,1);C=zeros(1,n+1);for p=1:n-1[Y,j]=max(abs(B(p:n,p)));C=B(p,:);B(p,:)=B(j+p-1,:);B(j+p-1,:)=C;for k=p+1:nm=B(k,p)/B(p,p);B(k,p:n+1)=B(k,p:n+1)-m*B(p,p:n+1);endendb=B(1:n,n+1);A=B(1:n,1:n);X(n)=b(n)/A(n,n);for q=n-1:-1:1X(q)=(b(q)-sum(A(q,q+1:n)*X(q+1:n)))/A(q,q);endelsedisp('请注意：因为RA=RB<n，所以方程组有无穷多解') endend输入：>> A=[10 -7 0 1;-3 2.099999 6 2;5 -1 5 -1;2 1 0 2];>> b=[8;5.900001;5;1];>> [RA,RB,n,X]=liezhu(A,b),H=det(A)输出：请注意：因为RA=RB，所以方程组有唯一解RA =4RB =4n =4X =0.0000-1.00001.00001.0000H =-762.0001第二题：建立列主元高斯消去法m文件（题一中已有）（1）输入：>> format compact>> A=[3.01 6.03 1.99;1.27 4.16 -1.23;0.987 -4.81 9.34];>> b=[1;1;1];>> [RA,RB,n,X]=liezhu(A,b),h=det(A),C=cond(A)输出：请注意：因为RA=RB，所以方程组有唯一解RA =3RB =n =3X =1.0e+03 *1.5926-0.6319-0.4936h =-0.0305C =3.0697e+04(2)输入：>> A=[3.00 6.03 1.99;1.27 4.16 -1.23;0.990 -4.81 9.34]; >> b=[1;1;1];>> [RA,RB,n,X]=liezhu(A,b),h=det(A)输出：请注意：因为RA=RB，所以方程组有唯一解RA =3RB =3n =3X =119.5273-47.1426-36.8403h =-0.4070第三题：输入：>> clear>> A=[10 7 8 7;7 5 6 5;8 6 10 9;7 5 9 10];>> b=[32 23 33 31]’;>> dA=det(A),lamda=eig(A),Ac2=cond(A,2)输出：dA =1.0000lamda =0.01020.84313.858130.28872.9841e+03下面分析误差性态：建立m文件：function Acp=pjwc(A,jA,b,jb,p)%Acp矩阵A的p条件数cond%pjwc：p范数解的误差性态分析%jA是A的近似矩阵jA=A+δA，jb=b+δbAcp=cond(A,p);dA=det(A);X=A\b;deltaA=jA-A;pndA=norm(deltaA,p);deltab=jb-b;pndb=norm(deltab,p);if pndb>0jX=A\jb;Pnb=norm(b,p);pnjx=norm(jX,p);deltaX=jX-X;pnjdX=norm(deltaX,p);jxX=pnjdX/pnjX;pnX=norm(X,p);xX=pnjdX/pnX;pndb=norm(deltab,p);xAb=pndb/pnb;pnbj=norm(jb,p);xAbj=pndb/pnbj;Xgxx=Acp*xAb;endif pndA>0jX=jA\b;deltaX=jX-X;pnX=norm(X,p);pnjdX=norm(deltaX,p);pnjX=norm(jX,p);jxX=pnjdX/pnjX;xX=pnjdX/pnX;pnjA=norm(jA,p);pnA=norm(A,p);pndA=norm(deltaA,p);xAbj=pndA/pnjA;xAb=pndA/pnA;Xgxx=Acp*xAb;endif (Acp>50)&(dA<0.1)disp('请注意：AX=b是病态的，A的p条件数Acp，A的行列式值dA，解X，近似解jX，解的相对误差xX，解的相对误差估计Xgxx，b或A的相对误差xAb依次如下：')Acp,dA,X',jX',xX',jxX',Xgxx',xAb',xAbj'elsedisp('请注意：AX=b是良态的，A的p条件数Acp，A的行列式值dA，解X，近似解jX，解的相对误差xX，解的相对误差估计Xgxx，b或A的相对误差xAb依次如下：')Acp,dA,X',jX',xX',jxX',Xgxx',xAb',xAbj'end输入：>> jA=[10 7 8.1 7.2;7.08 5.04 6 5;8 5.98 9.89 9;6.99 5 9 9.98];>> jb=b;p=2;>> Acp=pjwc(A,jA,b,jb,p)输出：请注意：AX=b是良态的，A的p条件数Acp，A的行列式值dA，解X，近似解jX，解的相对误差xX，解的相对误差估计Xgxx，b或A的相对误差xAb依次如下：Acp =2.9841e+03dA =1.0000ans =1.0000 1.0000 1.0000 1.0000 ans =-9.5863 18.3741 -3.2258 3.5240 xX =10.4661jxX =0.9842Xgxx =22.7396xAb =0.0076xAbj =0.0076Acp =2.9841e+03第四题：（1）输入：建立m文件:for n=2:6a=hilb(n);pnH(n-1)=cond(a,inf);endpnHn=2:6;plot(n,pnH);可见条件数随着n的增大而急剧增大（2）输入：>> n=2;H=hilb(n);>> x=(linspace(1,1,n))';>> b=H*x;>> [RA,RB,n,X]=gauss(H,b)输出：请注意：因为RA=RB，所以方程组有唯一解RA =2RB =2n =2X =1.00001.0000输入：>> r=b-H*X,deltax=X-x输出：r =deltax =1.0e-15 *0.4441-0.6661输入：>> n=3;H=hilb(n);>> x=(linspace(1,1,n))';>> b=H*x;>> [RA,RB,n,X]=gauss(H,b)>> r=b-H*X,deltax=X-x输出：X =1.00001.00001.0000r =1.0e-15 *0.2220deltax =1.0e-13 *-0.02000.1221-0.1255>> n=4;H=hilb(n);>> x=(linspace(1,1,n))'; >> b=H*x;>> [RA,RB,n,X]=gauss(H,b) >> r=b-H*X,deltax=X-xX =1.00001.00001.00001.0000r =1.0e-15 *-0.4441-0.1110deltax =1.0e-12 *-0.02220.2485-0.59800.3886>> n=5;H=hilb(n);>> x=(linspace(1,1,n))'; >> b=H*x;>> [RA,RB,n,X]=gauss(H,b) >> r=b-H*X,deltax=X-xX =1.00001.00001.00001.00001.0000r =1.0e-15 *0.22200.1110deltax =1.0e-11 *-0.00350.0524-0.19370.2591-0.1148>> n=6;H=hilb(n);>> x=(linspace(1,1,n))'; >> b=H*x;>> [RA,RB,n,X]=gauss(H,b) >> r=b-H*X,deltax=X-xX =1.00001.00001.00001.00001.0000r =1.0e-15 *0.22200.1110deltax =1.0e-11 *-0.00350.0524-0.19370.2591-0.1148>> n=7;H=hilb(n);>> x=(linspace(1,1,n))'; >> b=H*x;>> [RA,RB,n,X]=gauss(H,b) >> r=b-H*X,deltax=X-xX =1.00001.00001.00001.00001.00001.0000r =1.0e-15 *0.22200.1110deltax =1.0e-09 *-0.00080.0219-0.14820.3854-0.42540.1677>> n=8;H=hilb(n);>> x=(linspace(1,1,n))'; >> b=H*x;>> [RA,RB,n,X]=gauss(H,b) >> r=b-H*X,deltax=X-xX =1.00001.00001.00001.00001.00001.00001.00001.0000r =1.0e-15 *-0.2220-0.1110-0.1110deltax =1.0e-06 *-0.00000.0018-0.02360.1279-0.34420.4870-0.34660.0978>> n=9;H=hilb(n);>> x=(linspace(1,1,n))'; >> b=H*x;>> [RA,RB,n,X]=gauss(H,b) >> r=b-H*X,deltax=X-xX =1.00001.00001.00001.00001.00001.00001.00001.00001.0000r =1.0e-15 *0.44410.2220-0.22200.22200.2220-0.1110deltax =1.0e-04 *-0.00000.0002-0.00280.0197-0.07220.1471-0.16870.1017-0.0251>> n=10;H=hilb(n);>> x=(linspace(1,1,n))';>> b=H*x;>> [RA,RB,n,X]=gauss(H,b)>> r=b-H*X,deltax=X-xX =1.00001.00001.00001.00000.99991.00030.99961.00040.99981.0000r =1.0e-15 *0.4441-0.22200.2220-0.11100.11100.1110deltax =1.0e-03 *-0.00000.0001-0.00230.0205-0.09740.2669-0.43690.4214-0.22090.0485一、填空题(每空1分，共40分)1、中国共产党是_中国工人阶级___的先锋队，同时是中国人民和中华民族的先锋队，是_______中国特色社会主义事业___ 的领导核心。
数值分析课程实验设计——数值积分实习题
数值分析——数值积分实习题管理科学与工程学院 学号：1120140500 姓名：彭洋洋 一、计算实习题1.用不同数值方法计算积分：049xdx =-⎰.（1）取不同的步长h ，分别用复合梯形及复合辛普森求积计算积分，给出误差中关于h 的函数，并与积分精确值比较两个公式的精度，是否存在一个最小的h ，使得精度不能再被改善？ （2）用龙贝格求积计算完成问题（1） （3）用自适应辛普森积分，使其精度达到10-4解答：（1）取不同的步长，采用不同的公式，比较精度过程如下： 1.1 复合梯形公式及复合辛普森公式求解复合梯形公式：11*[()2()()]2n n k k hT f a f x f b -==++∑误差关于h 的函数：2(2)()**()12n a b R f h f ξ-=复合辛普森公式：111/201*[()4()2()()]6n n n k k k k hS f a f x f x f b --+===+++∑∑误差关于h 的函数：4(4)()*(/2)*()180n a bR f h f η-=1.2 复合梯形公式及复合辛普森公式Matlab 程序（2）用龙贝格求积计算完成问题（1） 2.1 龙贝格求积算法龙贝格求积公式也称为逐次分半加速法。
它是在梯形公式、辛普森公式和柯特斯公式之间的关系的基础上，构造出一种加速计算积分的方法。
作为一种外推算法，它在不增加计算量的前提下提高了误差的精度。
24133n n n S T T =- 21611515n n n C S S =- 26416363n n n R C C =-1221/201()22n n n k k h T T f x -+==+∑ ()(1)()11(4*)/(41)k m k k mm m m T T T +--=-- 1,2,...k = 2.2 龙贝格求积Matlab 程序画图程序设计①得到关于n各种公式求积的图表如下：对于梯形公式、辛普森公式n代表份数，龙贝格公式n表示从1开始的序列号②关于步长h 的各种公式求积的图表如下其中龙贝格序列步长()/2k h b a =-：观察两幅图表h 越小，精度越高。
数值分析第五章实习题答案
数值分析第五章实习题答案数值分析第五章实习题答案数值分析是一门研究如何使用计算机来解决数学问题的学科。
在数值分析的学习过程中，实习题是非常重要的一部分，通过实习题的练习，可以帮助我们巩固所学的知识，并且提高我们的解题能力。
本文将为大家提供数值分析第五章实习题的答案，希望对大家的学习有所帮助。
第一题：求下列方程的一个正根，并用二分法和牛顿法分别计算根的近似值。
方程：x^3 - 3x + 1 = 0解答：首先，我们可以通过绘制函数图像来初步估计方程的根的范围。
根据图像，我们可以大致确定根在区间[0, 2]之间。
接下来，我们使用二分法来计算根的近似值。
根据二分法的原理，我们将区间[0, 2]等分为两部分，然后判断根在哪一部分。
不断重复这个过程，直到找到根的近似值。
具体计算过程如下：- 将区间[0, 2]等分为两部分，得到中点x = 1。
- 计算方程在x = 1处的函数值f(1) = -1。
- 根据函数值的正负性，我们可以确定根在区间[1, 2]之间。
- 将区间[1, 2]等分为两部分，得到中点x = 1.5。
- 计算方程在x = 1.5处的函数值f(1.5) = 1.375。
- 根据函数值的正负性，我们可以确定根在区间[1, 1.5]之间。
- 重复以上步骤，直到找到根的近似值。
最终得到根的近似值为x ≈ 1.365。
接下来，我们使用牛顿法来计算根的近似值。
牛顿法是一种迭代法，通过不断逼近根的位置来计算根的近似值。
具体计算过程如下：- 选择初始近似值x0 = 1。
- 计算方程在x = 1处的函数值f(1) = -1。
- 计算方程在x = 1处的导数值f'(1) = 4。
- 利用牛顿法的迭代公式x1 = x0 - f(x0)/f'(x0)，我们可以得到x1 ≈ 1.333。
- 重复以上步骤，直到找到根的近似值。
最终得到根的近似值为x ≈ 1.365。
通过二分法和牛顿法，我们分别得到了方程x^3 - 3x + 1 = 0的一个正根的近似值为x ≈ 1.365。
数值分析计算实习作业一
数值分析计算实习题一学号：：院系：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] 上三次样条插值更精确。
数值分析计算实习题
插值法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]上三次样条插值更精确。
数值分析第五版计算实习题
弟二草插值法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].上.三次样条插值更精确，儿乎与原函数帀合。
数值分析计算实习题答案
数值分析计算实习题答案数值分析计算实习题答案数值分析是一门研究如何利用计算机对数学问题进行近似求解的学科。
在数值分析的学习过程中，实习题是一种重要的学习方式，通过实践来巩固理论知识，并培养解决实际问题的能力。
本文将为大家提供一些数值分析计算实习题的答案，希望能够帮助大家更好地理解和掌握数值分析的相关知识。
一、插值与拟合1. 已知一组数据点，要求通过这些数据点构造一个一次插值多项式，并求出在某一特定点的函数值。
答案：首先，我们可以根据给定的数据点构造一个一次插值多项式。
假设给定的数据点为(x0, y0), (x1, y1)，我们可以构造一个一次多项式p(x) = a0 + a1x，其中a0和a1为待定系数。
根据插值条件，我们有p(x0) = y0，p(x1) = y1。
将这两个条件代入多项式中，可以得到一个方程组，通过求解这个方程组，我们就可以确定a0和a1的值。
最后，将求得的多项式代入到某一特定点，就可以得到该点的函数值。
2. 已知一组数据点，要求通过这些数据点进行最小二乘拟合，并求出拟合曲线的表达式。
答案：最小二乘拟合是一种通过最小化误差平方和来找到最佳拟合曲线的方法。
假设给定的数据点为(x0, y0), (x1, y1)，我们可以构造一个拟合曲线的表达式y =a0 + a1x + a2x^2 + ... + anx^n，其中a0, a1, ..., an为待定系数。
根据最小二乘拟合原理，我们需要最小化误差平方和E = Σ(yi - f(xi))^2，其中yi为实际数据点的y值，f(xi)为拟合曲线在xi处的函数值。
通过求解这个最小化问题，我们就可以确定拟合曲线的表达式。
二、数值积分1. 已知一个函数的表达式，要求通过数值积分的方法计算函数在某一区间上的定积分值。
答案：数值积分是一种通过将定积分转化为数值求和来近似计算的方法。
假设给定的函数表达式为f(x)，我们可以将定积分∫[a, b]f(x)dx近似为Σwi * f(xi)，其中wi为权重系数，xi为待定节点。
实习题答案p50,p239(数值分析(第五版))
一1 求牛顿插值多项式、差商、插值及其误差估计的MATLAB主程序函数如下：function [y,R,A,C,L]=newdscg(X,Y,x,M)n=length(X); m=length(x);for t=1:mz=x(t); A=zeros(n,n);A(:,1)=Y';s=0.0; p=1.0; q1=1.0; c1=1.0;for j=2:nfor i=j:nA(i,j)=(A(i,j-1)- A(i-1,j-1))/(X(i)-X(i-j+1));endq1=abs(q1*(z-X(j-1)));c1=c1*j;endC=A(n,n);q1=abs(q1*(z-X(n)));for k=(n-1):-1:1C=conv(C,poly(X(k)));d=length(C);C(d)=C(d)+A(k,k);endy(k)= polyval(C, z);endR=M*q1/c1;L(k,:)=poly2sym(C);首先将名为newdscg.m的程序保存为M文件，然后在MA TLAB工作窗口输入程序>> syms M,X=[0.2,0.4,0.6,0.8,1.0];Y =[0.98,0.92,0.81,0.64,0.38];x=0.9;[y,R,A,C,P]=newdscg(X,Y,x,M)运行后输出插值y=f(0.9)及其误差限公式R，三阶牛顿插值多项式P及其系数向量C，差商的矩阵A如下:y =0.5239R =(7*M)/800000A =0.9800 0 0 0 00.9200 -0.3000 0 0 00.8100 -0.5500 -0.6250 0 00.6400 -0.8500 -0.7500 -0.2083 00.3800 -1.3000 -1.1250 -0.6250 -0.5208C =-0.5208 0.8333 -1.1042 0.1917 0.9800P =- (25*x^4)/48 + (5*x^3)/6 - (53*x^2)/48 + (23*x)/120 + 49/502:用拉格朗日插值法求其4次插值多项式的过程如下：先编写一个M文件来求其拉格朗日差值多项式，其程序如下：function f=Language(x,y,x0)syms t l;if(length(x)==length(y))n=length(x);elsedisp('x和y的维数不相等！');return;%检错endh=sym(0);for(i=1:n)l=sym(y(i));for(j=1:i-1)l=l*(t-x(j))/(x(i)-x(j));end;for(j=i+1:n)l=l*(t-x(j))/(x(i)-x(j));end;h=h+l;endsimplify(h);if(nargin==3)f=subs(h,'t',x0);%计算插值点的函数值?elsef=collect(h);f=vpa(f,6);%将插值多项式的系数化成6位精度的小数?end然后将Language.m保存为M文件,然后在MA TLAB工作窗口输入程序:>> x=[0.2 0.4 0.6 0.8 1.0];>> y=[0.98 0.92 0.81 0.64 0.38];>> f=Language(x,y,0.9)其运行结果为：f =0.5239然后在输入如下程序可得其4次差值多项式>> f=Language(x,y)其运行结果为：f =- 0.520833*t^4 + 0.833333*t^3 - 1.10417*t^2 + 0.191667*t + 0.98最后输入>> plot(x,y)可得如下图二1.牛顿迭代法：dfun.m:function fun = dfun( x )%求原方程的导数f′(x)%UNTITLED4 Summary of this function goes here% Detailed explanation goes heresyms tf=diff(fun0(t));fun=subs(f,t,x);endfunction fun = fun0(x)%UNTITLED3 Summary of this function goes here% Detailed explanation goes herefun=x^2-3*x+2-exp(x);%所要求的原方程endnewton1.mfunction [y ,k ] = newton1 (x0 ,n ,derta)%y 为迭%代序列,k 为迭代次数,x0 初始位置,n 为最大迭代次数,derta 为%容许误差%UNTITLED5 Summary of this function goes here% Detailed explanation goes herek = 1 ;y = [x0] ;t = x0 - fun0(x0)./dfun(x0) ;while abs (t - x0) >=dertay=[y ,t ] ;x0 = t ;k = k +1 ;t = x0 - fun0(x0)./dfun(x0) ;if (k -1) > nerror ('n is full') ;endEnd1.其结果：>> [y ,k ] = newton1(1,40,0.000000001)y =1.0000 0.2689 0.2575 0.2575k =42.不动点迭代：程序如下：f=inline('(x^2+2-exp(x))/3');x=feval(f,0.5);disp(x);Eps=1E-8;i=1;while 1i=i+1;x=feval(f,x);disp(x);if abs(x-x0)<Epsdisp(i)break;endend运行结果如下：0.2004262430999570.2727490650983750.2536071565841300.2585503762649360.2572656363350940.2575989851621900.2575124545148320.2575349136152510.25752908416795693.斯特芬森迭代法程序format long;f=inline('x-((x^2+2-exp(x))/3-x)^2/((((x^2+2-exp(x))/3)^2+2-exp((x^2+2-exp(x))/3))/ 3-2*(x^2+2-exp(x))/3+x)');disp('x=');x=feval(f,0.5);disp(x);Eps=1E-8;i=1;while 1x0=x;i=i+1;x=feval(f,x);disp(x);if abs(x-x0)<Epsdisp(i)break;endend运行结果0.2586844275657910.2575303177198080.2575302854398610.257530285439861 4。
数值分析计算实习第一题
直接用定义: ������������(������������)2 = ‖������������‖2‖������������−1‖2
求 A 的条件数很繁琐,需要先进行化简：
首先：
由于 A 是对称矩阵，
‖������������‖2 = �������������max(������������������������������������)
说明 ：
1. 在所用的算法中，凡是要给出精度水平的ε，都取 ������������=10−12。
2. 选择算法的时候应使矩阵 A 的所有零元素都不存储。
3. 打印以下内容：
(1)算法设计方案和思路。
(2)全部源程序。
(3)特征值������������1,������������501,������������������������,������������������������������������(������������=1,2,⋯,39)以及������������������������������������������������(������������)2, det������������的值(采用 e 型输出实型数，并 至少显示 12 位有效数字)。
λi[16] -2.533970311130E+00 λi[38] 8.648666065193E+00
λi[17] -2.003230769563E+00 λi[39] 9.254200344575E+00
λi[18] -1.503557611227E+00 cond(A)2 1.925204273903E+03
λi[19] -9.935586060080E-01 det(A) 2.772786141752E+118
数值分析(第五版)计算实习题第四章作业
第四章：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轴')。
p178-179 实习题答案(数值分析(第五版))
1.列主元法（178页实习1题）函数liezhu(A,b).m的头文件的程序如下：n=length(b);detA=det(A);d=0;x=zeros(n,1);c=zeros(1,n);if abs(det(A))<=epserror('系数矩阵是奇异的'); %所求方程组系数矩阵必须是是非奇异的！return;end%按列选主元for i=1:n-1max=abs(A(i,i));m=i;for j=i+1:nif max<abs(A(j,i))max=abs(A(j,i));m=j;endend%换行if m~=ifor k=i:nc(k)=A(i,k);A(i,k)=A(m,k);A(m,k)=c(k);endd=b(i);b(i)=b(m);b(m)=d;end%消元for k=i+1:nfor j=i+1:nA(k,j)=A(k,j)-A(i,j)*A(k,i)/A(i,i);endb(k)=b(k)-b(i)*A(k,i)/A(i,i);A(k,i)=0;endEnd%回代求解x(n)=b(n)/A(n,n);for i=n-1:(-1):1sum=0;for j=i+1:nsum=sum+A(i,j)*x(j);endx(i)=(b(i)-sum)/A(i,i);endend其结果如下：>> A=[10,-7,0,1;-3,2.099999,6,2;5,-1,5,-1;2,1,0,2];>> b=[8;5.900001;5;1];>> [x,detA]=liezhu(A,b)x =0.0000-1.00001.00001.0000detA =-762.00012.LU分解法：（178页实习1题）LUFJ.m文件中函数LUFJ(A,b)的程序如下：function LUFJ(A,b) %A为系数矩阵，b为右端项矩阵%UNTITLED2 Summary of this function goes here% Detailed explanation goes here[m,n]=size(A); %初始化矩阵A，b，L和Un=length(b);L=eye(n,n);U=zeros(n,n);U(1,1:n)=A(1,1:n); %开始进行LU分解L(2:n,1)=A(2:n,1)/U(1,1);for k=2:nU(k,k:n)=A(k,k:n)-L(k,1:k-1)*U(1:k-1,k:n);L(k+1:n,k)=(A(k+1:n,k)-L(k+1:n,1:k-1)*U(1:k-1,k))/U(k,k); endL %输出L矩阵U %输出U矩阵y=zeros(n,1); %开始解方程组Ux=yy(1)=b(1);for k=2:ny(k)=b(k)-L(k,1:k-1)*y(1:k-1);endx=zeros(n,1);x(n)=y(n)/U(n,n);for k=n-1:-1:1x(k)=(y(k)-U(k,k+1:n)*x(k+1:n))/U(k,k); endfor k=1:nfprintf('x[%d]=%f\n',k,x(k));endend运行结果如下：>> A=[10,-7,0,1;-3,2.099999,6,2;5,-1,5,-1;2,1,0,2];>> b=[8;5.900001;5;1];>> LUFJ(A,b)L =1.0e+006 *0.0000 0 0 0-0.0000 0.0000 0 00.0000 -2.5000 0.0000 00.0000 -2.4000 0.0000 0.0000U =1.0e+007 *0.0000 -0.0000 0 0.00000 -0.0000 0.0000 0.00000 0 1.5000 0.57500 0 0 0.0000x[1]=-0.000000x[2]=-1.000000x[3]=1.000000x[4]=1.0000003.雅可比迭代法：（179页实习3题）Jacobi.m文件中函数Jacobi(A,b,eps)的程序如下：function Jacobi(A,b,eps) %A为系数矩阵，b为后端项矩阵，epe为精度[m,n]=size(A);D=diag(diag(A)); %求矩阵DL=D-tril(A); %求矩阵LU=D-triu(A); %求矩阵Utemp=1;x=zeros(m,1);k=0;while abs(max(x)-temp)>epstemp=max(abs(x));k=k+1; %记录循环次数x=-inv(D)*(L+U)*x+inv(D)*b; %雅克比迭代公式endfor k=1:nfprintf('x[%d]=%f\n',k,x(k));end运行结果如下：>> A=[10 -1 2 0;0 8 -1 3;2 -1 10 0;-1 3 -1 11];>> b=[-11;-11;6;25];>> Jacobi(A,b,0.00005)x[1]=-1.466538x[2]=-2.358692x[3]=0.657438x[4]=2.8424524 .Gauss-Seidel迭代程序（179页实习3题）Gauss-Seidel.m文件中函数Gauss-Seidel(A,b,eps)的程序如下：function Gauss_Seidel(A,b,eps) %A为系数矩阵，b为后端项矩阵，epe为精度[m,n]=size(A);D=diag(diag(A)); %求矩阵DL=D-tril(A); %求矩阵LU=D-triu(A); %求矩阵Utemp=1;x=zeros(m,1);k=0;while abs(max(x)-temp)>epstemp=max(abs(x));k=k+1; %记录循环次数x=inv(D-L)*U*x+inv(D-L)*b; %Gauss_Seidel的迭代公式endfor k=1:nfprintf('x[%d]=%f\n',k,x(k));end其运行结果如下：>> A=[10 -1 2 0;0 8 -1 3;2 -1 10 0;-1 3 -1 11];>> b=[-11;-11;6;25];>> Gauss_Seidel(A,b,0.005) x[1]=-1.466538x[2]=-2.358692x[3]=0.657438x[4]=2.842452。
数值分析计算实习题目一
数值分析计算实习题目一SY0905308 徐捷设有501501⨯的矩阵123499500501..........a b c b a b c c b a b c A cb a bc cb a b cba ⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦其中()()()0.11.640.024sin 0.20.641,2,,501;0.16,0.064.ii a i i ei b c =--=⋅⋅⋅==-矩阵A 的特征值()1,2,,501i i λ=⋅⋅⋅满足12501,λλλ<<⋅⋅⋅<1501||min ||s i i λλ≤≤=, 试求:1.1501,,s λλλ的值2. A 的与数5011140k kλλμλ-=+最接近的特征值(1,2,,39)k i k λ=3. A 的(谱范数)条件数2()cond A 和行列式det A .1. 算法的设计方案本题的核心算法是幂法、带原点平移的幂法、反幂法和LU 分解法，要点在于选择算法时，应使A 的所有零元素都不存储。
故算法设计的思路如下，第一步，对A 使用幂法(Powermethod)，可得A 的按模最大的特征值，记为a λ； 第二步，对A 使用带有原点平移的幂法，令平移量1a p λ=，可得另一端点的特征值记为b λ；第三步，比较a λ与b λ的大小，根据条件12501,λλλ<<⋅⋅⋅<可知{}1min ,a b λλλ=，{}501max ,a b λλλ=；第四步，对A 使用反幂法(Inversepowermethod)，可得A 的按模最小的特征值s λ（使用LU 杜立特尔分解法）第五步，根据5011140k kλλμλ-=+计算出k μ,然后利用带有原点平移的反幂法求得k i λ,其中平移量k k p μ=,反幂法运算39次,可得2239,,,i i i λλλ ;第六步, 根据定义,非奇异的实对称矩阵A 的谱范数条件数()12||ncond A λλ=,其中1n λλ和分别是矩阵A 的模为最大和模为最小的特征值，对于本题，则有{}15012|max ||,|||()||s cond A λλλ=；第七步，由LU 分解可知，A LU =，可得5011det det()det()det()(1501)n iii A LU L U u i ====⋅=≤≤∏。
数值分析第五版计算实习题
数值分析计算实习题第二章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: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)*df(i);enddisp('4次牛顿插值多项式');P4=vpa(collect((sum(d))),5) %P4即为4次牛顿插值多项式,并保留小数点后5位数pp=csape(x1,y1, 'variational');%调用三次样条函数q=pp.coefs;disp('三次样条函数');for i=1:4S=q(i,:)*[(x-x1(i))^3;(x-x1(i))^2;(x-x1(i));1];S=vpa(collect(S),5)endx2=0.2:0.08:1.08;dot=[1 2 11 12];figureezplot(P4,[0.2,1.08]);hold ony2=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三次样条函数x∈[0.2,0.4]时， S = - 1.3393*x^3 + 0.80357*x^2 - 0.40714*x + 1.04 x∈[0.4,0.6]时，S = 0.44643*x^3 - 1.3393*x^2 + 0.45*x + 0.92571x∈[0.6,0.8]时，S = - 1.6964*x^3 + 2.5179*x^2 - 1.8643*x + 1.3886 x∈[0.8,1.0]时，S =2.5893*x^3 - 7.7679*x^2 + 6.3643*x - 0.80571输出图如下2-3〔1〕程序：clear;clc;x1=[0 1 4 9 16 25 36 49 64];y1=[0 1 2 3 4 5 6 7 8];%插值点n=length(y1);a=ones(n,2);a(:,2)=-x1';c=1;for i=1:nc=conv(c,a(i,:));endq=zeros(n,n);r=zeros(n,n+1);for i=1:n[q(i,:),r(i,:)]=deconv(c,a(i,:));%wn+1/(x-xk) end 三次样条插值曲线4次牛顿插值曲线Dw=zeros(1,n);for i=1:nDw(i)=y1(i)/polyval(q(i,:),x1(i));%系数endp=Dw*q;syms x L8;for i=1:nL8(i)=p(n-i+1)*x^(i-1);enddisp('8次拉格朗日插值');L8=vpa(collect((sum(L8))),5)xi=0:64;yi=polyval(p,xi);figureplot(xi,yi,x1,y1,'r*');hold ontitle('8次拉格朗日插值');结果如下：8次拉格朗日插值L8 =- 3.2806e-10*x^8 + 6.7127e-8*x^7 - 5.4292e-6*x^6 + 0.00022297*x^5 - 0.0049807*x^4 + 0.060429*x^3 - 0.38141*x^2 + 1.3257*x输出图如下：第五章4-1〔3〕程序：clc;clear;y= (x) sqrt(x).*log(x);a=0;b=1;tol=1e-4;p=quad(y,a,b,tol);fprintf('采用自适应辛普森积分结果为： %d \n', p); 结果如下：采用自适应辛普森积分结果为： -4.439756e-01第九章9-1(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由结果可知，步长越小，结果越精确。
数值分析计算实习题
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]上三次样条插值更精确。