数值分析计算实习题
数值分析计算实习题二
《数值分析》计算实习题二算法设计方案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，该函数可对系数矩阵相同// 而右端向量不同的多个方程组同时求解。
数值分析上机实习题
2019-2020 第1学期数值分析上机实习题总目标：会算，要有优化意识。
（以下程序要求以附件1例题代码格式给出）1. 对给定的线性方程组Ax b =进行迭代求解。
（1）给出Jacobi 迭代的通用程序。
（2）给出Gauss-Seidel 迭代的通用程序。
调用条件：系数矩阵A ，右端项b ，初值0x ，精度要求ε。
输出结果：方程组的近似解。
给定线性方程组211122241125x --⎛⎫⎛⎫ ⎪ ⎪= ⎪ ⎪ ⎪ ⎪---⎝⎭⎝⎭，和122711122215x -⎛⎫⎛⎫ ⎪ ⎪= ⎪ ⎪ ⎪ ⎪⎝⎭⎝⎭，取初值0x 为0， 分别利用Jacobi 迭代和G-S 迭代进行求解，观察并解释其中的数学现象。
2. 利用紧凑格式（即直接分解法或逐框运算法）对给定的矩阵A 进行Doolittle 分解，并用其求线性方程组的解。
调用条件：矩阵A 。
输出结果：单位下三角矩阵L 和上三角矩阵U 。
给定矩阵1112A ⎛⎫= ⎪⎝⎭，利用以下算法：1）将A 作Doolittle 分解11A LU =，2）令211A U L =，并对2A 作Doolittle 分解222A L U =，3）重复2）的过程令11n n n A U L --=，并对n A 作Doolittle 分解n n n A L U =，2,3,4,n =， 观察n L ，n U ，n A 的变化趋势，思考其中的数学现象。
3. 给定函数21(),12511f x x x -≤+≤=，取164,8,n =，用等距节点21,i i n x =-+ 0,1,,1i n =+对原函数进行多项式插值和五次多项式拟合，试画出插值和拟合曲线，并给出数学解释。
4. 给出迭代法求非线性方程()0f x =的根的程序。
调用条件：迭代函数()x ϕ，初值0x输出结果：根的近似值k x 和迭代次数k给定方程32()10f x x x =--=，用迭代格式1k x +=0 1.5x =附近的根，要使计算结果具有四位有效数字，利用估计式*1||1||k k k L x x x x L -≤---，或估计式*10||1||kk L x x x x L-≤--来判断需要的迭代次数，分别需要迭代多少次？两者是否有冲突？5. 利用数值求积算法计算()ba f x dx ⎰。
数值分析(第五版)计算实习题第五章作业教学资料
>> 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，所以方程组有唯一解
ans =
-9.5863 18.3741 -3.2258 3.5240
xX =
10.4661
jxX =
0.9842
Xgxx =
22.7396
xAb =
0.0076
xAbj =
0.0076
Acp =
2.9841e+03
第四题：
（1）输入：
建立m文件:
forn=2:6
a=hilb(n);
pnH(n-1)=cond(a,inf);
RA =
3
RB =
3
n =
3
X =
1.0e+03 *
1.5926
-0.6319
-0.4936
h =
-0.0305
C =
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)
>> r=b-H*X,deltax=X-x
输出：
X =
数值分析第五章实习题答案
数值分析第五章实习题答案数值分析第五章实习题答案数值分析是一门研究如何使用计算机来解决数学问题的学科。
在数值分析的学习过程中，实习题是非常重要的一部分，通过实习题的练习，可以帮助我们巩固所学的知识，并且提高我们的解题能力。
本文将为大家提供数值分析第五章实习题的答案，希望对大家的学习有所帮助。
第一题：求下列方程的一个正根，并用二分法和牛顿法分别计算根的近似值。
方程：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。
数值分析计算实习题
插值法1.下列数据点的插值x 0 1 4 9 16 25 36 49 64y 0 1 2 3 4 5 6 7 8可以得到平方根函数的近似，在区间[0,64] 上作图.(1)用这9 个点作8 次多项式插值Ls(x).(2)用三次样条( 第一边界条件)程序求S(x).从得到结果看在[0,64] 上，哪个插值更精确；在区间[0,1] 上，两种插值哪个更精确解：(1) 拉格朗日插值多项式，求解程序如下syms x l;x1=[0 1 4 9 16 25 36 49 64]; y1=[0 1 2 3 4 5 6 7 8]; n=length(x1);Ls=sym(0);for i=1:nl=sym(y1(i));for k=1:i-1l=l*(x-x1(k))/(x1(i)-x1(k));endfor k=i+1:n l=l*(x-x1(k))/(x1(i)-x1(k));endLs=Ls+l; endLs=simplify(Ls) % 为所求插值多项式Ls(x).输出结果为Ls = -/*xA2+95549/72072*x-1/00*xA8-2168879/0*xA4+19/0*xA7+657859/*xA3+33983/ 0*xA5-13003/00*xA6(2) 三次样条插值，程序如下x1=[0 1 4 9 16 25 36 49 64];y1=[0 1 2 3 4 5 6 7 8];x2=[0:1:64];y3=s plin e(x1,y1,x2);p=po Iyfit(x2,y3,3); % 得到三次样条拟合函数S=p(1)+p(2)*x+ p(3)*x^2+p(4)*xA3 % 得到S(x) 输出结果为：S =/6464-2399/88*x+/1984*xA2+2656867/624*xA3⑶ 在区间［0,64］上，分别对这两种插值和标准函数作图,Plot(x2,sqrt(x2),'b',x2,y2,'r',x2,y3,'y')蓝色曲线为y="X函数曲线，红色曲线为拉格朗日插值函数曲线，黄色曲线为三次样条插值曲线可以看到蓝色曲线与黄色曲线几乎重合，因此在区间[0,64] 上三次样条插值更精确。
数值分析第五版计算实习题
弟二草插值法3.卜列数据点的插值可以得到平方根函数的近似，在区间064］上作图。
(1〉用这9个点做8次多项式插值Q x)。
(2)用三次样条(第一边界条件)程岸求S(X)。
从得到结果石在［0.64］ 1：・哪个插值更粘确：在区间［0,1］ I：•两种插值哪个更精确?(1) 8次多项式插值：(1)8次多项式插值:首先建立新的M-file:输入如卜代码(此为拉格朗口插值的功能函数)并保存function f=Language(x,y,x0)%求Li知数据点的拉格朗Fl插值多项式%己知数据点的x坐标向量：x%已知数据点的y坐标向量:y%插值的x坐标：x0%求得的拉格朗H插值多项式或在X0处的插值：fsyms t;ifi(lcngth(x)=length(y))n=length(x);elsedisp(*x和y的维数不相等!)；return;end %检错tbr(i=l:n)i=y(i);fbr(j=1:i-l)l=l*(t-x(j))/(x(i)-x(j))；end;for(j=i-M:n)end;for(j=i+l:n) l=l*(t-x(j))/(x(i)-x(j))； end;simplify(f);if(i==n) if|nargin=3)f=subs(C't\xO);else f=collcct(f);f=vpa(f,6);endendend再建立新的M-file：输入：clear;x=[0 1 49 16 25 36 49 64];y=[0：l：8]；%计算拉格朗口基丞数%计算拉格朗ri插值函数%化简%计算插值点的曲数值%将插值多项式展开%将插值多项式的系数化成6位精度的小数f=Uinguage(x,y) 运行得到f=1.32574*1-381410*t A2+.604294e-1 *t A3+.222972e-3 *t A5-.542921 e-5*t A6+.671268e・7T7・.328063e・9T8・.498071 e-2*t A4 这就是8次多项式插值L s(x)= 1.32574怜.381410*t A2+.604294e-1 *t A3+.222972e-3 *t A5-.542921 e-5*t A6+.671268e-7*t A7-.328063e-9*t A8-.498071 e-2*t A4. (2)三次样条插值：建立新的M-filc：输入：clear;x=[0 I 49 1625 36 4964];尸[0:8];t=[0:0.1:64];Y=t.A(0.5);O=Language(x,y)f= 1,32574*t-.381410*t.A2+.604294e-1 *t.A3+.222972e-3*t.A5-.542921 e・5*(. W+.671268e-7*t.A7-.328063e-9*t.A8-.498071 e-2 *t.A4;S=interp l(x,y,t.'spline,);plol(x,y,o；(・YY.lf.'b'」S'g:');grid;运行程序得到如下图：从结果屮很明显可以看出在[0.64].上.三次样条插值更精确，儿乎与原函数帀合。
北航数值分析-实习作业1(C语言详细注释)
《数值分析》计算实习作业《一》北航第一题 设有501501⨯的矩阵⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡=501500499321a bc b a b cc b a b ccb a bc c b a b c b a A其中.064.0,16.0);501,2,1(64.0)2.0sin()024.064.1(1.0-===--=c b i e i i a i i 矩阵的特征值)501,,2,1( =i i λ满足||min ||,501150121i i s λλλλλ≤≤=<<<试求1. 5011，λλ和s λ的值2. 的与数4015011λλκλμ-+=k 最接近的特征值)39,,2,1( =K κλi3. 的（谱范数）条件数2)A (cond 和行列式A det 要求1. 算法的设计方案（A 的所有零元素都不能存储）2. 全部源程序（详细注释）。
变量为double ，精度-1210=ε，输出为e 型12位有效数字3. 特征值s 5011，，λλλ和)39,,2,1( =K κλi 以及A cond det ,)A (2的值4. 讨论迭代初始向量的选取对计算结果的影响，并说明原因解答：1. 算法设计对于s λ满足||min ||5011i i s λλ≤≤=，所以s λ是按模最小的特征值，直接运用反幂法可求得。
对于5011，λλ，一个是最大的特征值，一个是最小的特征值，不能确定两者的绝对值是否相等，因此必须首先假设||||5011λλ≠，然后运用幂法，看能否求得一个特征值，如果可以求得一个，证明A 是收敛的，求得的结果是正确的，然后对A 进行带原点平移的幂法，偏移量是前面求得的特征值，可以求得另一个特征值，最后比较这两个特征值，较大的特征值是501λ，较小的特征值就是1λ。
如果在假设的前提下，无法运用幂法求得按模最大的特征值，即此时A 不收敛，则需要将A 进行带原点平移的幂法，平移量可以选取1，再重复上述步骤即可求得两个特征值。
北航数值分析计算实习题目一 幂法 反幂法 求矩阵特征值
《数值分析》计算实习题目第一题：1. 算法设计方案（1）1λ，501λ和s λ的值。
1)首先通过幂法求出按模最大的特征值λt1，然后根据λt1进行原点平移求出另一特征值λt2，比较两值大小，数值小的为所求最小特征值λ1，数值大的为是所求最大特征值λ501。
2)使用反幂法求λs ，其中需要解线性方程组。
因为A 为带状线性方程组，此处采用LU 分解法解带状方程组。
（2）与140k λλμλ-5011=+k 最接近的特征值λik 。
通过带有原点平移的反幂法求出与数k μ最接近的特征值 λik 。
（3）2cond(A)和det A 。
1）1=nλλ2cond(A)，其中1λ和n λ分别是按模最大和最小特征值。
2）利用步骤（1）中分解矩阵A 得出的LU 矩阵，L 为单位下三角阵,U 为上三角阵，其中U 矩阵的主对角线元素之积即为det A 。
由于A 的元素零元素较多，为节省储存量，将A 的元素存为6×501的数组中，程序中采用get_an_element()函数来从小数组中取出A 中的元素。
2.全部源程序#include <stdio.h>#include <math.h>void init_a();//初始化Adouble get_an_element(int,int);//取A 中的元素函数double powermethod(double);//原点平移的幂法double inversepowermethod(double);//原点平移的反幂法int presolve(double);//三角LU 分解int solve(double [],double []);//解方程组int max(int,int);int min(int,int);double (*u)[502]=new double[502][502];//上三角U 数组double (*l)[502]=new double[502][502];//单位下三角L 数组double a[6][502];//矩阵Aint main(){int i,k;double lambdat1,lambdat2,lambda1,lambda501,lambdas,mu[40],det;double lambda[40];init_a();//初始化Alambdat1=powermethod(0);lambdat2=powermethod(lambdat1);lambda1=lambdat1<lambdat2?lambdat1:lambdat2;lambda501=lambdat1>lambdat2?lambdat1:lambdat2;presolve(0);lambdas=inversepowermethod(0);det=1;for(i=1;i<=501;i++)det=det*u[i][i];for (k=1;k<=39;k++){mu[k]=lambda1+k*(lambda501-lambda1)/40;presolve(mu[k]);lambda[k]=inversepowermethod(mu[k]);}printf("------------所有特征值如下------------\n");printf("λ=%1.11e λ=%1.11e\n",lambda1,lambda501);printf("λs=%1.11e\n",lambdas);printf("cond(A)=%1.11e\n",fabs(lambdat1/lambdas));printf("detA=%1.11e \n",det);for (k=1;k<=39;k++){printf("λi%d=%1.11e ",k,lambda[k]);if(k % 3==0) printf("\n");} delete []u;delete []l;//释放堆内存return 0;}void init_a()//初始化A{int i;for (i=3;i<=501;i++) a[1][i]=a[5][502-i]=-0.064;for (i=2;i<=501;i++) a[2][i]=a[4][502-i]=0.16;for (i=1;i<=501;i++) a[3][i]=(1.64-0.024*i)*sin(0.2*i)-0.64*exp(0.1/i); }double get_an_element(int i,int j)//从A中节省存储量的提取元素方法{if (fabs(i-j)<=2) return a[i-j+3][j];else return 0;}double powermethod(double offset)//幂法{int i,x1;double u[502],y[502];double beta=0,prebeta=-1000,yita=0;for (i=1;i<=501;i++)u[i]=1,y[i]=0;//设置初始向量u[]for (int k=1;k<=10000;k++){yita=0;for (i=1;i<=501;i++) yita=sqrt(yita*yita+u[i]*u[i]);for (i=1;i<=501;i++) y[i]=u[i]/yita;for (x1=1;x1<=501;x1++){u[x1]=0;for (int x2=1;x2<=501;x2++)u[x1]=u[x1]+((x1==x2)?(get_an_element(x1,x2)-offset):get_an_element(x1,x2))*y[x2] ;}prebeta=beta;beta=0;for (i=1;i<=501;i++) beta=beta+ y[i]*u[i];if (fabs((prebeta-beta)/beta)<=1e-12) {printf("offset=%f lambda=%f err=%e k=%d\n",offset,(beta+offset),fabs((prebeta-beta)/beta),k);break;};//输出中间过程,包括偏移量,误差,迭代次数}return (beta+offset);}double inversepowermethod(double offset)//反幂法{int i;double u[502],y[502];double beta=0,prebeta=0,yita=0;for (i=1;i<=501;i++)u[i]=1,y[i]=0; //设置初始向量u[]for (int k=1;k<=10000;k++){yita=0;for (i=1;i<=501;i++) yita=sqrt(yita*yita+u[i]*u[i]);for (i=1;i<=501;i++) y[i]=u[i]/yita;solve(u,y);prebeta=beta;beta=0;for (i=1;i<=501;i++) beta=beta+ y[i]*u[i];beta=1/beta;if (fabs((prebeta-beta)/beta)<=1e-12) {printf("offset=%f lambda=%f err=%e k=%d\n",offset,(beta+offset),fabs((prebeta-beta)/beta),k);break;};//输出中间过程,包括偏移量,误差,迭代次数}return (beta+offset);int presolve(double offset)//三角LU分解{int i,k,j,t;double sum;for (k=1;k<=501;k++)for (j=1;j<=501;j++){u[k][j]=l[k][j]=0;if (k==j) l[k][j]=1;} //初始化LU矩阵for (k=1;k<=501;k++){for (j=k;j<=min(k+2,501);j++){sum=0;for (t=max(1,max(k-2,j-2)) ; t<=(k-1) ; t++)sum=sum+l[k][t]*u[t][j];u[k][j]=((k==j)?(get_an_element(k,j)-offset):get_an_element(k,j))-sum;}if (k==501) continue;for (i=k+1;i<=min(k+2,501);i++){sum=0;for (t=max(1,max(i-2,k-2));t<=(k-1);t++)sum=sum+l[i][t]*u[t][k];l[i][k]=(((i==k)?(get_an_element(i,k)-offset):get_an_element(i,k))-sum)/u[k][k];}}return 0;}int solve(double x[],double b[])//解方程组{int i,t;double y[502];double sum;y[1]=b[1];for (i=2;i<=501;i++){sum=0;for (t=max(1,i-2);t<=i-1;t++)sum=sum+l[i][t]*y[t];y[i]=b[i]-sum;}x[501]=y[501]/u[501][501];for (i=500;i>=1;i--){sum=0;for (t=i+1;t<=min(i+2,501);t++)sum=sum+u[i][t]*x[t];x[i]=(y[i]-sum)/u[i][i];}return 0;}int max(int x,int y){return (x>y?x:y);}int min(int x,int y){return (x<y?x:y);}3.计算结果结果如下图所示:部分中间结果：给出了偏移量(offset)，误差(err)，迭代次数(k)4.讨论迭代初始向量的选取对计算结果的影响,并说明原因使用u[i]=1（i=1,2,...,501）作为初始向量进行迭代，可得出以上结果。
数值分析计算实习二
void atra(double a[N][N], int n);//矩阵转置
void aadd(double a[N], double b[N], double c[N], int n, int co);//向量加法
void aaadd(double a[N][N], double b[N][N], double c[N][N], int n, int co);//矩阵加法
结果：
1、 拟上三角矩阵 如下。
2、迭代最终矩阵如下，迭代次数为14次。
3、矩阵所有特征值如下。
从上至下10个特征值，其中左列为实部，右列为虚部，可知有两对复特征值。
4、对应于实特征值的特征向量如下
问题发现与讨论
1、双步位移QR迭代过程中迭代矩阵 始终为拟上三角矩阵。
2、双步位移QR方法与基本QR方法比较
amul(wr, ur, D, m+1);
aaadd(A, D, A, m+1, -1);
amul(ur, pr, D, m+1);
aaadd(A, D, A, m+ห้องสมุดไป่ตู้, -1);
}
}
}
k++;
if (k > 10000)//算法失败判定
{
printf("算法失败\n");
break;
}
}
for (int i = 0; i < N; i++)
doubleaneiji(double a[N], double b[N], int n);//向量内积
void aamul(double a[N], double b[N][N], double c[N], int n);//矩阵向量乘法
数值分析计算实习第一题
直接用定义: ������������(������������)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。
数值分析实验题目
实验一：函数插值与数据拟合1.1实验目的（1）由函数()f x 的1n +个节点处函数值得出n 次Lagrange 插值函数； （2）由函数()f x 的1n +个节点处函数值得出n 次Newton 插值函数；（3）由函数()f x 的个1n +节点处函数值得出Hermite 插值函数或分段三次Hermite 函数； （4）由未知函数的离散数据(),1,2,,i f x i n = 得出最小二乘拟合函数。
1.2 实验原理（1）00,()()nnj n k k j j kk jx x L x y x x ==≠-=⋅-∑∏（2）00100120101011()()[,]()[,,]()()[,,,]()()()n n n N x f x f x x x x f x x x x x x x f x x x x x x x x x -=+-+--++---（3）222100,0,0,1[12()]()n n nn i i n k k k kk i i k i i k i i k k i k i k i x x x x H y x x y x x x x x x x x +==≠=≠=≠⎧⎫⎛⎫⎛⎫--⎪⎪'=⋅--+⋅-⎨⎬ ⎪ ⎪---⎝⎭⎝⎭⎪⎪⎩⎭∑∑∏∏（4）()()()()()()()()()()()()0001000101111101001122,,,,,,,,,,,,n n n n n n n n n na y a y a y S a a a a φφφφφφφφφφφφφφφφφφφφφφφφφ⎛⎫⎛⎫⎛⎫⎪ ⎪ ⎪⎪⎪ ⎪= ⎪ ⎪ ⎪⎪ ⎪⎪ ⎪ ⎪⎝⎭⎝⎭⎝⎭=++++1.3实验内容1.3,===用Lagrange23,===用Newton3. 给定函数21(),551f x x x=-≤≤+，取插值节点5,0,1,2,,10k x k k =-+= ，求分段三次Hermite 插值函数，并在一个坐标系中画出两函数的图形。
数值分析试题及答案
数值分析试题及答案一、单项选择题（每题3分，共30分）1. 线性代数中，矩阵A的逆矩阵记作（）。
A. A^TB. A^-1C. A^+D. A*答案：B2. 插值法中，拉格朗日插值多项式的基函数是（）。
A. 多项式B. 指数函数C. 正弦函数D. 余弦函数答案：A3. 在数值积分中，梯形规则的误差是（）阶的。
A. O(h^2)B. O(h^3)C. O(h)D. O(1/h)答案：A4. 求解线性方程组时，高斯消元法的基本操作不包括（）。
A. 行交换B. 行乘以非零常数C. 行加行D. 行除以非零常数答案：D5. 非线性方程f(x)=0的根的迭代法中，收敛的必要条件是（）。
A. f'(x)≠0B. f'(x)=0C. |f'(x)|<1D. |f'(x)|>1答案：C6. 利用牛顿法求解非线性方程的根时，需要计算（）。
A. 函数值B. 函数值和导数值C. 函数值和二阶导数值D. 函数值、一阶导数值和二阶导数值答案：B7. 矩阵的特征值和特征向量是（）问题中的重要概念。
A. 线性方程组B. 特征值问题C. 线性规划D. 非线性方程组答案：B8. 在数值分析中，条件数是衡量矩阵（）的量。
A. 稳定性B. 可逆性C. 正交性D. 稀疏性答案：A9. 利用龙格现象说明，高阶插值多项式在区间端点附近可能产生（）。
A. 振荡B. 收敛C. 稳定D. 单调答案：A10. 雅可比迭代法和高斯-塞德尔迭代法都是求解线性方程组的（）方法。
A. 直接B. 迭代C. 精确D. 近似答案：B二、填空题（每题4分，共20分）11. 线性代数中，矩阵A的行列式记作________。
答案：det(A) 或 |A|12. 插值法中，牛顿插值多项式的基函数是________。
答案：差商13. 在数值积分中，辛普森规则的误差是________阶的。
答案：O(h^4)14. 求解线性方程组时，迭代法的基本思想是从一个初始近似解出发，通过不断________来逼近精确解。
数值分析计算实习题
1.下列数据点的插值x 0 1 4 9 16 25 36 49 64y 0 1 2 3 4 5 6 7 8可以得到平方根函数的近似，在区间[0,64]上作图.(1)用这9个点作8次多项式插值Ls(x).(2)用三次样条(第一边界条件)程序求S(x).从得到结果看在[0,64]上，哪个插值更精确；在区间[0,1]上，两种插值哪个更精确？解：(1)拉格朗日插值多项式，求解程序如下syms x l;x1=[0 1 4 9 16 25 36 49 64];y1=[0 1 2 3 4 5 6 7 8];n=length(x1);Ls=sym(0);for i=1:nl=sym(y1(i));for k=1:i-1l=l*(x-x1(k))/(x1(i)-x1(k));endfor k=i+1:nl=l*(x-x1(k))/(x1(i)-x1(k));endLs=Ls+l;endLs=simplify(Ls) %为所求插值多项式Ls(x).输出结果为Ls =-/*x^2+95549/72072*x-1/00*x^8-2168879/0*x^4+19/0*x^7+657859/*x^3+33983/ 0*x^5-13003/00*x^6(2)三次样条插值，程序如下x1=[0 1 4 9 16 25 36 49 64];y1=[0 1 2 3 4 5 6 7 8];x2=[0:1:64];y3=spline(x1,y1,x2);p=polyfit(x2,y3,3); %得到三次样条拟合函数S=p(1)+p(2)*x+p(3)*x^2+p(4)*x^3 %得到S(x)输出结果为：S =/6464-2399/88*x+/1984*x^2+2656867/624*x^3(3)在区间[0,64]上，分别对这两种插值和标准函数作图，plot(x2,sqrt(x2),'b',x2,y2,'r',x2,y3,'y')蓝色曲线为y=函数曲线,红色曲线为拉格朗日插值函数曲线，黄色曲线为三次样条插值曲线可以看到蓝色曲线与黄色曲线几乎重合，因此在区间[0,64]上三次样条插值更精确。
数值分析计算实习题（二）
数值分析计算实习题（二）数值分析计算实习题（二）SY1004114 全昌彪一：算法设计方案概述：本题采用fortran90语言编写程序，依据题目要求，采用带双步位移QR分解法求出所给矩阵的所有特征值，并求出相应于其实特征值的特征向量，以及相关需要给出的中间结果。
1、矩阵的A的初始化（赋值）：利用子函数initial（a,n）来实现，返回n×n 维二维数组a。
2、A矩阵的拟上三角化：利用子函数hessenberg(a,n)，在对矩阵进行QR分解前进行拟上三角化，这样可以提高计算效率，减少计算量，返回A矩阵的相似矩阵Hessenberg阵A(n-1)。
3、对A(n-1)进行带双步位移QR分解得出Ｃｍ及A矩阵的所有特征值，这一步利用了两个子函数eigenvalue(a,n,lamda,lamdai)和qrresolve(b,c,m)带双步位移QR分解可以加速收敛。
每次QR分解前先进行判断，若可以直接得到矩阵的特征值，则对矩阵直接降阶处理；若不可以，则进行QR分解，这样就进一步减少了计算量，提高了计算效率。
考虑到矩阵A可能有复特征值，采用两个一维数组lamda(n)及lamdai(n)分别存储其实部和虚部。
在双步位移处理及降阶过程中，被分解的矩阵Ａk(m ×m)及中间矩阵M k(m×m)的维数随m不断减少而降阶，于是引入了动态矩阵C(m×m)和B(m×m)分别存储，在使用前，先声明分配内存，使用结束后立即释放内存。
返回A(n-1)经双步位移QR分解后的矩阵及Ａ矩阵的所有特征值。
4、特征向量的求解：采用子函数eigenvector(a,lamda)实现求解A矩阵的属于实特征值的特征向量。
核心算法为高斯列主元消去法，(A-λI)x=b,b=0,回代过程令x(10)=1，即可求出对应于每一实特征值的特征向量的各个元素。
5、相关输出结果：所有数据均采用e型输出，数据保留到小数点后12位。
数值分析实习第一题
{
break;
}
}
}
}
void VECT_EVA(double Vect_new[], double Vect_bef[])//向量赋值
{
int i;
for(i = 0; i < N; i++)
{
Vect_new[i] = Vect_bef[i];
}
}
double RANK_MAX(double M_Array[], unsigned char n)//找最大值
{
if(i >= N)
{
break;
}
Sum_Ele[0] = 0;
Sum_Ele[1] = i - R;
Sum_Ele[2] = k - S;
for(t = RANK_MAX(Sum_Ele,3); t <= k - 1; t++)
{
LU_Mat[M(i,k)][k] -= LU_Mat[M(i,t)][t]*LU_Mat[M(t,k)][k];
}
LU_Mat[M(i,k)][k] = LU_Mat[M(i,k)][k]/LU_Mat[S][k];
}
}
}
}
void SOLVE_EQU(double Matrix_A[][N], double Targ_b[], double Ret_Un[])//解方程
{
int i, j, t;
double Sol_Ele[2];
//In_Beta = VECT_T_VEC(In_Vect_y, In_Vect_u);
In_Beta = 1.0/VECT_T_VEC(In_Vect_y, In_Vect_u);
数值分析第九章计算实习题答案昆工
数值分析第九章计算实习题答案昆工(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、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
《数值分析》计算实习题姓名：学号：班级：第二章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-1)*df(i);endP4=vpa(sum(d),5) %P4即为4次牛顿插值多项式,并保留小数点后5位数pp=csape(x1,y1, 'variational');%调用三次样条函数q=pp.coefs;q1=q(1,:)*[(x-.2)^3;(x-.2)^2;(x-.2);1];q1=vpa(collect(q1),5)q2=q(1,:)*[(x-.4)^3;(x-.4)^2;(x-.4);1];q2=vpa(collect(q2),5)q3=q(1,:)*[(x-.6)^3;(x-.6)^2;(x-.6);1];q3=vpa(collect(q3),5)q4=q(1,:)*[(x-.8)^3;(x-.8)^2;(x-.8);1];q4=vpa(collect(q4),5)%求解并化简多项式2、运行结果P4 =0.98*x - 0.3*(x - 0.2)*(x - 0.4) - 0.625*(x - 0.2)*(x - 0.4)*(x - 0.6) - 0.20833*(x - 0.2)*(x - 0.4)*(x - 0.8)*(x - 0.6) + 0.784q1 =- 1.3393*x^3 + 0.80357*x^2 - 0.40714*x + 1.04q2 =- 1.3393*x^3 + 1.6071*x^2 - 0.88929*x + 1.1643 q3 =- 1.3393*x^3 + 2.4107*x^2 - 1.6929*x + 1.4171 q4 =- 1.3393*x^3 + 3.2143*x^2 - 2.8179*x + 1.86293、问题结果4次牛顿差值多项式4()P x = 0.98*x - 0.3*(x - 0.2)*(x - 0.4) - 0.625*(x - 0.2)*(x- 0.4)*(x - 0.6) - 0.20833*(x - 0.2)*(x - 0.4)*(x - 0.8)*(x - 0.6) + 0.784三次样条差值多项式()Q x0.10.20.30.40.50.60.70.80.910.40.50.60.70.80.911.1323232321.33930.803570.40714 1.04,[0.2,0.4]1.3393 1.60710.88929 1.1643,[0.4,0.6]1.3393 2.4107 1.6929 1.4171,[0.6,0.8]1.3393 3.21432.8179 1.8629,[0.8,1.0]x x x x x x x x x x x x x x x x ⎧-+-+∈⎪-+-+∈⎪⎨-+-+∈⎪⎪-+-+∈⎩第三章1、程序代码Clear;clc;x=[0 0.1 0.2 0.3 0.5 0.8 1];y=[1 0.41 0.5 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,'c--',x,y1,'r:',x,y2,'y-.')p3=polyfit(x,y,2)%观察图像，类似抛物线，故用二次多项式拟合。
y3=polyval(p3,x);plot(x,y,'c--',x,y1,'r:',x,y2,'y-.',x,y3,'k--')%画出四种拟合曲线2、运行结果p1 =-6.6221 12.8147 -4.6591 0.9266 p2 =2.8853 -12.3348 16.2747 -5.2987 0.9427 p3 =3.1316 -1.2400 0.73563、问题结果三次多项式拟合P1=32-6.622112.8147 4.65910.9266x x x +-+四次多项式拟合P2=4322.885312.334816.2747 5.29870.9427x x x x -+-+ 二次多项式拟合P3=23.1316 1.24000.7356x x -+第四章1、程序代码1)建立函数文件f.m: function y=fun(x); y=sqrt(x)*log(x); 2)编写程序：a. 利用复化梯形公式及复化辛普森公式求解：Clear;clc;h=0.001;%h 为步长,可分别令h=1,0.1,0.01,0.001 n=1/h;t=0;s1=0;s2=0; for i=1:n-1 t=t+f(i*h); endT=h/2*(0+2*t+f(1));T=vpa(T,7) %梯形公式0.10.20.30.40.50.60.70.80.9100.511.522.53for i=0:n-1s1=s1+f(h/2+i*h);endfor i=1:n-1s2=s2+f(i*h);endS=h/6*(0+4*s1+2*s2+f(1));S=vpa(S,7) %辛普森公式a’复化梯形公式和复化辛普生公式程序代码也可以是：Clear;clc;x=0:0.001:1; %h为步长,可分别令h=1,0.1,0.01,0.001y=sqrt(x).*log(x+eps);T=trapz(x,y);T=vpa(T,7)（只是h=1的运行结果不一样，T=1.110223*10^(-16)，而其余情况结果都相同）Clear;clc;f=inline('sqrt(x).*log(x)',x);S=quadl(f,0,1);S=vpa(S,7)b.利用龙贝格公式求解：Clear;clc;m=14;%m+1即为二分次数h=2;for i=1:mh=h/2;n=1/h;t=0;for j=1:n-1t=t+f(j*h);endT(i)=h/2*(0+2*t+f(1));%梯形公式endfor i=1:m-1for j=m:i+1T(j)=4^i/(4^i-1)*T(j)-1/(4^i-1)*T(j-1);%通过不断的迭代求得T（j)，即T表的对角线上的元素。
endendvpa(T(m),7)2、运行结果T =-0.4443875S =-0.4444345ans =-0.44444143、问题结果b. 利用龙贝格公式求解：通过15次二分，得到结果：-0.4444414第五章1、程序代码（1）LU分解解线性方程组：Clear;clc;A=[10 -7 0 1-3 2.099999 6 25 -1 5 -12 1 0 2];b=[8;5.900001;5;1];[m,n]=size(A);L=eye(n);U=zeros(n);flag='ok';for i=1:nU(1,i)=A(1,i);endfor r=2:nL(r,1)=A(r,1)/U(1,1);endfor i=2:nfor j=i:nz=0;for r=1:i-1z=z+L(i,r)*U(r,j);endU(i,j)=A(i,j)-z;endif abs(U(i,i))<epsflag='failure'return;endfor k=i+1:nm=0;for q=1:i-1m=m+L(k,q)*U(q,i);endL(k,i)=(A(k,i)-m)/U(i,i);endendLUy=L\b;x=U\ydetA=det(L*U)（2）列主元消去法：function x = gauss(A,b);A=[10 -7 0 1;-3 2.099999 6 2;5 -1 5 -1;2 1 0 2];b=[8;5.900001;5;1];[n,n] = size(A);x = zeros(n,1);Aug = [A,b]; %增广矩阵for k = 1:n-1[piv,r] = max(abs(Aug(k:n,k))); %找列主元所在子矩阵的行r r = r + k - 1; % 列主元所在大矩阵的行if r>ktemp=Aug(k,:);Aug(k,:)=Aug(r,:);Aug(r,:)=temp;endif Aug(k,k)==0, error(‘对角元出现0’), end% 把增广矩阵消元成为上三角for p = k+1:nAug(p,:)=Aug(p,:)-Aug(k,:)*Aug(p,k)/Aug(k,k); endend% 解上三角方程组A = Aug(:,1:n); b = Aug(:,n+1);x(n) = b(n)/A(n,n);for k = n-1:-1:1x(k)=b(k);for p=n:-1:k+1x(k) = x(k)-A(k,p)*x(p);endx(k)=x(k)/A(k,k);enddetA=det(A)2、运行结果1）LU分解解线性方程组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 =-0.0000-1.00001.00001.0000detA =-762.00012）列主元消去法detA =762.0001ans =0.0000-1.00001.00001.00003、问题结果1）LU分解解线性方程组L=10003100 0.5250000010 0.224000000.961⎛⎫ ⎪-⎪ ⎪-⎪-⎝⎭U=10701006 2.3 00150000055749998.499 000 5.079998908 -⎛⎫ ⎪ ⎪ ⎪ ⎪⎝⎭x=（0.0000，−1.0000，1.0000，1.0000）T detA=-762.0012)列主元消去法x=（0.0000，−1.0000，1.0000，1.0000）T detA=762.001第六章1、程序代码（1）Jacobi迭代Clear;clc;n = 6; %也可取n=8,10H = hilb(n);b = H * ones(n, 1);e = 0.00001;for i = 1:nif H(i, i)==0 '对角元为零，不能求解' returnendendx = zeros(n, 1);k = 0;kend = 10000;r = 1;while k<=kend & r>ex0 = x;for i = 1:ns = 0;for j = 1:i - 1s = s + H(i, j) * x0(j);endfor j = i + 1:ns = s + H(i, j) * x0(j);endx(i) = b(i) / H(i, i) - s / H(i, i); endr = norm(x - x0, inf);k = k + 1;endif k>kend '迭代不收敛，失败'else '求解成功'xkend（2）SOR迭代1）程序代码function s = SOR(n, w);H = hilb(n);b = H*ones(n, 1);e = 0.00001;for i = 1:nif H(i,i)==0 ‘对角线为零，不能求解’returnendendx = zeros(n, 1);k = 0;kend = 10000;r = 1;while k<=kend & r>ex0 = x;for i = 1:ns = 0;for j = 1:i - 1s = s + H(i, j) * x(j);endfor j = i + 1:ns = s + H(i, j) * x0(j);endx(i) = (1 - w) * x0(i) + w / H(i, i) * (b(i) - s); endr = norm(x - x0, inf);k = k + 1;endif k>kend '迭代不收敛，失败'else '求解成功'xend2）从命令窗口中分别输入：SOR（6,1）SOR（8,1）SOR（10,1）SOR（6,1.5）SOR（8,1.5）SOR（10，1.5）2、运行结果Jacobi迭代：ans =迭代不收敛，失败SOR迭代：第七章 1、程序代码 (1)不动点迭代法1）建立函数文件：g.m function f=g(x)f(1)=20/(x^2+2*x+10); 2)建立函数文件：bdd.mfunction [y, n] = bdd(x, eps) if nargin==1 eps=1.0e-8; elseif nargin<1 error return endx1 = g(x); n = 1;while (norm(x1-x)>=eps)&&(n<=10000) x = x1;x1 = g(x);n = n + 1;endy = x;n3)从命令窗口输入：bdd(0)（2）牛顿迭代clear;clc;format long;m=8; %m为迭代次数，可分别令m=2,4,6,8,10x=sym('x');f=sym('x^3+2*x^2+10*x-20');df=diff(f,x);FX=x-f/df; %牛顿迭代公式Fx=inline(FX);disp('x=');x1=0.5;disp(x1);Eps=1E-8;k=0;while 1x0=x1;k=k+1;x1=feval(Fx,x1); %将x1代入牛顿迭代公式替代x1 disp(x1); %在屏幕上显示x1if k==mbreak;endendk,x12、运行结果（1）不动点迭代法>> bdd(0)n =25ans =1.3688（2）牛顿迭代x=21.4666666666666671.3715120138059211.3688102226338951.3688081078226671.3688081078213731.3688081078213731.368808107821373k =8x1 =1.3688081078213733、问题结果(1)不动点迭代法x=1.3688 n=25 收敛太慢（2）牛顿迭代初值取0迭代次数k=8时，Welcome To Download !!!欢迎您的下载，资料仅供参考！。