数值分析课程实验设计——数值积分实习题
数值分析课程实验设计——数值线性代数实习题
数值分析课程实验设计——数值线性代数
实习题
1. 实验目的
本实验的主要目的是进一步加深对数值线性代数的理解，熟悉
常见矩阵分解方法，并在此基础上解决实际问题。
2. 实验内容
本次实验将任务分为两个部分，分别是矩阵分解与求解线性方
程组。
2.1 矩阵分解
首先，我们需要熟悉三种常见的矩阵分解：QR分解、LU分解
和奇异值分解。
我们需要通过Python语言实现这三种分解方法，
并利用这些方法解决实际问题。
2.2 求解线性方程组
其次，我们需要学会用矩阵分解的方法来求解线性方程组。
我
们将通过两个例子来进行说明，并利用Python语言实现这些方法。
3. 实验要求
本次实验要求熟悉矩阵分解的基本方法，在此基础上解决实际问题；能够运用多种方法来求解线性方程组，并分析比较它们的优缺点。
4. 实验总结
本次实验通过矩阵分解和求解线性方程组两个部分的学习，巩固了我们对于数值线性代数的知识，并在实际问题的解决中得到了应用。
感谢老师的指导，我们会在今后的学习中持续探索数值分析方面的知识。
数值分析上机实习题
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 ⎰。
数值分析计算实习题
数值分析计算实习题-CAL-FENGHAI-(2020YEAR-YICAI)_JINGBIAN《数值分析》计算实习题姓名：学号：班级：第二章1、程序代码Clear;clc;x1=[ ];y1=[ ];n=length(y1);c=y1(:);for j=2:n %求差商for i=n:-1:jc(i)=(c(i)-c(i-1))/(x1(i)-x1(i-j+1));endendsyms x df d;df(1)=1;d(1)=y1(1);for i=2:n %求牛顿差值多项式df(i)=df(i-1)*(x-x1(i-1));d(i)=c(i-1)*df(i);endP4=vpa(sum(d),5) %P4即为4次牛顿插值多项式,并保留小数点后5位数pp=csape(x1,y1, 'variational');%调用三次样条函数q=;q1=q(1,:)*[^3;^2;;1];q1=vpa(collect(q1),5)q2=q(1,:)*[^3;^2;;1];q2=vpa(collect(q2),5)q3=q(1,:)*[^3;^2;;1];q3=vpa(collect(q3),5)q4=q(1,:)*[^3;^2;;1];q4=vpa(collect(q4),5)%求解并化简多项式2、运行结果P4 =*x - *(x - *(x - - *(x - *(x - *(x - - *(x - *(x - *(x - *(x - + q1 =- *x^3 + *x^2 - *x +q2 =- *x^3 + *x^2 - *x + q3 =- *x^3 + *x^2 - *x + q4 =- *x^3 + *x^2 - *x +3、问题结果4次牛顿差值多项式4()P x = *x - *(x - *(x - - *(x - *(x - *(x - - *(x - *(x - *(x - *(x - +三次样条差值多项式()Q x0.10.20.30.40.50.60.70.80.910.40.50.60.70.80.911.1323232321.33930.803570.40714 1.04,[0.2,0.4]1.3393 1.60710.88929 1.1643,[0.4,0.6]1.3393 2.4107 1.6929 1.4171,[0.6,0.8]1.3393 3.21432.8179 1.8629,[0.8,1.0]x x x x x x x x x x x x x x x x ⎧-+-+∈⎪-+-+∈⎪⎨-+-+∈⎪⎪-+-+∈⎩第三章1、程序代码Clear;clc; x=[0 1]; y=[1 ];p1=polyfit(x,y,3)%三次多项式拟合 p2=polyfit(x,y,4)%四次多项式拟合 y1=polyval(p1,x);y2=polyval(p2,x);%多项式求值plot(x,y,'c--',x,y1,'r:',x,y2,'y-.')p3=polyfit(x,y,2)%观察图像，类似抛物线，故用二次多项式拟合。
数值分析实验题目及解答
内容包括：实验题目1：算法的数值稳定性实验实验题目2：LU分解实验实验题目3：三次样条插值外推样条实验实验题目4：第二类Fredholm 积分方程实验实验题目5：M级显式R_K法实验题目：算法的数值稳定性实验实验内容：计算积分()10()d 1515nx I n x a x==+⎰ (n=1,2,…,20) 易得到下面递推公式()()11I n aI n n=--+并有估计式()()()()11111I n a n a n <<+++计算方法：算法一:采用下面递推公式计算:()()11I n aI n n =--+()1,2,,20n = 取初值()1160ln ln 15a I a +==算法二: 采用下面递推公式计算:()()111I n I n a n ⎡⎤-=-+⎢⎥⎣⎦()20,19,,1n =结果分析：（分析哪个好哪个不好，原因是什么） 我觉得算法二比较好, 原因一:根据式()()()()11111I n a n a n <<+++得知,I(n)不可能小于零,而算法一的计算结果有部分结果小于零。
原因二:对算法一记初始误差ε0=/I 0-I(0)/>0;则εn =/I n -I(n)/=a/I n-1-I(n-1)/=a n*ε0由此可知,当n=20时, ε20把ε0放大了a 20倍，其结果造成严重的。
而对于算法二^^11n n a εε-=,…, ^^01n n aεε=,尽管有初始误差^20ε,但随着计算的进程,这个误差的影响不断减小。
附：源程序：（把源程序附上） 算法一程序: >> format long>> a=15;I=log(16/15); for n=1:20 nI=-a*I+1/n end算法二程序: >> format long>> a=15;I=31/10080; >> for n=20:-1:1 n II=1/a*(-I+1/n); End。
数值分析实验二 数值积分1
数值分析实验二 数值积分组号 班级 学号 姓名 分数一：实验目的1、掌握用复化Simpson 公式，复化梯形求积公式计算积分的方法。
2、掌握用龙贝格Romberg 积分公式计算积分的方法。
3、掌握用高斯-勒让德Gauss-Legendre 公式计算积分法。
4、通过实例了解三种方法的联系与区别,并会利用适当的方法计算某函数在某个区间的积分值。
二：实验内容及基本知识介绍（1）复化Simpson 求积公式的原理：将区间[a,b]分为n 等分，在每个子区间[1,k k x x +]上采用辛普森公式：()()462b a a b s f a f f b -⎡+⎤⎛⎫=++ ⎪⎢⎥⎝⎭⎣⎦，若记122k k x h x +=+，则得： ()()()()()111110246k kn bx ax k n k k k n k I f x dx f x dxh f x f x f x R f +-=-++===⎡⎤⎛⎫=+++⎢⎥ ⎪⎢⎥⎝⎭⎣⎦∑⎰⎰∑。
则可以记得：()()()()()⎥⎥⎦⎤⎢⎢⎣⎡++⎪⎪⎭⎫ ⎝⎛+=⎥⎥⎦⎤⎢⎢⎣⎡+⎪⎪⎭⎫ ⎝⎛+=∑∑∑-=-=+-=++1010211012124646n k n k k k n k k k k n b f x f x f a f h x f x f x f h S，称为复化辛普森求积公式。
其余项为：()()()()4141,,1802n n n kkkk k h h R f I S f x x ηη-+=⎛⎫=-=-∈ ⎪⎝⎭∑，此外，由于nS中求积系数均为正数，故知复化辛普森公式计算稳定。
（2）复化梯形求积公式的原理：将区间[a,b]分为n 等分，分点,,2,1,0,,n k nab h kh a x k =-=+=在每个子区间[]()1,,1,0,1-=+n k x x k k 上采用梯形公式()()[]b f a f a b T +-=2计算，则得：()()()()[]()f R x f x f h dx x f dx x f I n n k k k n k x x bak k++===∑∑⎰⎰-=+-=+11121，记：()()[]()()()⎥⎦⎤⎢⎣⎡++=+=∑∑-=-=+`11101222n k k n k k k n b f x f a f h x f x f h T ，称为复化梯形公式。
数值分析计算实习题答案
数值分析计算实习题答案数值分析计算实习题答案数值分析是一门研究如何利用计算机对数学问题进行近似求解的学科。
在数值分析的学习过程中，实习题是一种重要的学习方式，通过实践来巩固理论知识，并培养解决实际问题的能力。
本文将为大家提供一些数值分析计算实习题的答案，希望能够帮助大家更好地理解和掌握数值分析的相关知识。
一、插值与拟合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为待定节点。
数值分析上机实验题参考
数值分析论文数值积分 一、问题提出选用复合梯形公式，复合Simpson 公式，Romberg 算法，计算I = dx x ⎰-4102sin 4 ()5343916.1≈II =dx x x ⎰1sin ()9460831.0,1)0(≈=I fI = dx xe x⎰+1024 I =()dx x x ⎰++1211ln 二、要求编制数值积分算法的程序；分别用两种算法计算同一个积分，并比较其结果；分别取不同步长()/ a b h -=n ，试比较计算结果（如n = 10, 20等）； ﹡给定精度要求ε，试用变步长算法，确定最佳步长﹡。
三、目的和意义深刻认识数值积分法的意义； 明确数值积分精度与步长的关系；根据定积分的计算方法，可以考虑二重积分的计算问题引言一、数值求积的基本思想实际问题当中常常需要计算积分，有些数值方法。
如微分方程和积分方程的求解，也都和积分计算相联系。
依据人们熟悉的微积分基本原理，对于积分I=⎰a b f(x)dx,只要找到被积函数f(x)和原函数F(x),便有下列牛顿-莱布尼茨公式：I=⎰a b f(x)dx=F(b)-F(a).但实际使用这种求积方法往往有困难，因为大量的被积函数，诸如x xsin，2xe-等，其原函数不能用初等函数表达，故不能用上述公式计算。
另外，当f(x)是由测量或数值计算给出的一张数据表时，牛顿-莱布尼茨公式也不能直接运用，因此有必要研究积分的数值计算问题。
二、数值积分代数精度数值求积方法是近似方法，为要保证精度，我们自然希望求积公式能对“尽可能多”的函数准确成立，就提出了所谓代数精度的概念。
如果某个求积公式对次数不超过m的多项式均能准确成立，但对m+1次多项式就不能准确成立，则称该求积公式具有m次代数精度。
三、复合求积公式为了提高精度，通常可以把积分区间分成若干子区间（通常是等分），再在每个子区间用低阶求积公式，即复化求积法，比如复化梯形公式与复化辛普森公式。
数值分析实验课——数值积分与数值微分实验
数值分析实验报告四数值积分与数值微分实验（2学时）一 实验目的1．掌握复化的梯形公式、辛扑生公式等牛顿-柯特斯公式计算积分。
2. 会用高斯公式计算积分。
3. 掌握数值微分的计算方法。
二 实验内容1．分别用复化的梯形公式和辛扑生公式计算积分。
⎰90dx x M=4 2．给定下列表格值利用四点式（n=3）求)50()50('''f f 和的值。
三 实验步骤（算法）与结果1复化的梯形公式：()()()()01121222bM M a h h h f x dx f f f f f f b a h M -=+++++-=⎰基于上述公式，进行编程，程序代码如下：#include"stdio.h"#include"math.h"#define M 4int main(){float a=9,b=0;int i;float h;float x[M+1];float f=0;h=(a-b)/M;for(i=0;i<M+1;i++){x[i]=b+h*i;}for(i=0;i<M;i++){f=f+h/2*((sqrt(x[i]))+(sqrt(x[i+1])));}printf("\n%f\n",f);return 0;}运行结果如下：复化的辛扑生公式：()()()()01351246243323bM M a M h h f x dx f f f f f f h f f f f b ah M--=+++++++++++-=⎰基于上述公式，进行编程计算，程序代码如下： #include"stdio.h"#include"math.h"#define M 4int main(){float a=9,b=0;int i;float h;float x[M+1];float f=0;h=(a-b)/M;for(i=0;i<M+1;i++){x[i]=b+h*i;}for(i=0;i<M/2;i++){f=f+(h/3)*((sqrt(x[2*i]))+4*(sqrt(x[2*i+1]))+(sqrt(x[2*i+2]))); }printf("\n%f\n",f);return 0;}运算结果如下：2一阶和二阶数值微分公式四点式：()()()()'00123''00123211118926112302466f x y y y y h f x y y y y h ≈-+-+≈-+-基于上述公式，进行编程计算，程序代码如下： #include"stdio.h"int main(){float x[4]={50,55,60,65};float y[4]={1.6990,1.7404,1.7782,1.8129};float f1,f2;float h=x[1]-x[0];f1=(1/(6*h))*(-11*y[0]+18*y[1]-9*y[2]+2*y[3]); f2=(1/(6*h*h))*(12*y[0]-30*y[1]+24*y[2]-6*y[3]);printf("f1=%f",f1);printf("\nf2=%f",f2);return 0;}运算结果如下：四实验收获与教师评语。
数值分析计算实习题
插值法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]上三次样条插值更精确。
数值分析上机实习题
数值分析上机实习题第2章插值法1. 已知函数在下列各点的值为试⽤四次⽜顿插值多项式）（x p 4及三次样条韩式)(S x (⾃然边界条件)对数据进⾏插值。
⽤图给出(){}10,11,1,0,08.02.0,,x i =+=i x y i i ，）（x p 4及)(x S Python 代码import numpy as npimport matplotlib.pyplot as pltfrom matplotlib.font_manager import FontPropertiesfont_set = FontProperties(fname=r"c:\windows\fonts\simsun.ttc",size=12) #求⽜顿n 次均差 def qiujuncha(x,f,n): for i in range(1,n): for j in range(4,i-1,-1):f[j]= (f[j] - f[j-1])/(x[j]-x[j-i]) #根据⽜顿多项式求值 def niudun(x,f,x1): sum = f[0]; tmp = 1;for i in range(1,5): tmp *= (x1-x[i-1]) sum = sum + f[i]*tmp return sum#⽜顿插值画图 def drawPic(x,f):x1 = np.linspace(0.2, 1, 100) plt.plot(x1, niudun(x,f,x1))plt.title(u"⽜顿四次插值",fontproperties=font_set) plt.xlabel(u"x 轴",fontproperties=font_set) plt.ylabel(u"y 轴", fontproperties=font_set) plt.show() def qiu_h(x,h): n = len(x) -1 for i in range(n): print(i)h[i] = x[i+1]-x[i]#⾃然边界条件下的三次样条插值求Mdef qiu_m(h,f,o,u,d):n = len(h)o[0] = 0u[n] = 0d[n] = d[0] = 0a = []for i in range(1,n):u[i] = h[i-1]/(h[i-1]+h[i])for i in range(1,n):o[i] = h[i]/(h[i-1]+h[i])for i in range(1,n-1):d[i] = 6*(f[i+1]-f[i])/(h[i-1]+h[i])t = [0 for i in range(5)]t[0] =2t[1] = o[0]a.append(t)for i in range(1,n):t = [0 for i in range(5)]t[i - 1] = u [i + 1]t[i] = 2t[i + 1] = o [i + 1]a.append(t)t = [0 for i in range(5)]t[n - 1] = u[n]t[n] = 2a.append(t)tmp = np.linalg.solve(np.array(a),np.array(d))m = []for i in range(5):m.append(tmp[i])return m#根据三次条插值函数求值def yangtiao(x1,m,x,y,h,j):returnm[j]*(x[j+1]-x1)**3/(6*h[j])+m[j+1]*(x1-x[j])**3/(6*h[j])+(y[j]-m[j]*h[j]**2/6)*(x[j+1]-x1)/h[j] +(y[j+1]-m[j+1]*h[j]**2/6)*(x1-x[j])/h[j] def main():x = [0.2, 0.4, 0.6, 0.8, 1.0]y = [0.98, 0.92, 0.81, 0.64, 0.38]f = y[:]f1 = y[:]h = [0.2,0.2,0.2,0.2]u = [0 for n in range(5)]d = [0 for n in range(5)]o = [0 for n in range(5)] qiujuncha(x,f,4) qiujuncha(x,f1,2)m = qiu_m(h,f1,o,u,d) x1 = np.linspace(0.2, 0.4, 10)p1= plt.plot(x1, yangtiao(x1,m,x,y,h,0),color='red') x1 = np.linspace(0.4, 0.6, 10)plt.plot(x1, yangtiao(x1, m, x, y, h, 1),color='red') x1 = np.linspace(0.6, 0.8, 10)plt.plot(x1, yangtiao(x1, m, x, y, h, 2),color='red') x1 = np.linspace(0.8, 1.0, 10)plt.plot(x1, yangtiao(x1, m, x, y, h, 3),color='red') x1 = np.linspace(0.2, 1.0, 40)p2 = plt.plot(x1,niudun(x,f,x1),color='green') plt.xlabel(u"x 轴", fontproperties=font_set) plt.ylabel(u"y 轴",fontproperties=font_set) plt.title("三次样条插值和⽜顿插值")plt.legend(labels=[u'三次样条插值',u'⽜顿插值'],prop=font_set,loc="best") plt.show() main()实验结果运⾏结果可得插值函数图（如图1-1），4次⽜顿插值函数）（x p 4和三次样条插值函数)(x S 如下：)6.0(*)4.0(*)2.0(625.0)4.0(*)2.0(*3.098.0)(4-------=x x x x x x x P 98.0)8.0(*)6.0(*)4.0(*)2.0(*20833.0+-----x x x x]4.0,2.0[),2.0(467.4)4.0(9.4)2.0(167.1)(S 3∈-+-+-=x x x x x]6.0,4.0[),4.0(113.4)6.0(6467.4)4.0(575.1)6.0(167.1)(S 33∈-+-+----=x x x x x x ]8.0,6.0[),6.0(2.3)8.0(113.4)6.0(575.1)(S 3∈-+-+--=x x x x x]0.1,8.0[),8.0(9.1)0.1(2.3)(S ∈-+-=x x x x图1-1三次样条插值和⽜顿插值图2.在区间[-1,1]上分别取n = 10,20⽤两组等距节点对龙格函数做多项式插值三次样条插值，对每个n值画出插值函数及图形。
数值分析实验(4)
页脚内容1实验四 数值积分与数值微分专业班级：信计131班 姓名：段雨博 学号：2013014907 一、实验目的1、熟悉matlab 编程。
2、学习数值积分程序设计算法。
3、通过上机进一步领悟用复合梯形、复合辛普森公式，以及用龙贝格求积方法计算积分的原理。
二、实验题目 P1371、用不同数值方法计算积分049xdx =-⎰。
（1）取不同的步长h .分别用复合梯形及复合辛普森求积计算积分，给出误差中关于h 的函数，并与积分精确值比较两个公式的精度，是否存在一个最小的h ，使得精度不能再被改善？（2）用龙贝格求积计算完成问题（1）。
三、实验原理与理论基础1.1复合梯形公式及其复合辛普森求解[]()()()11101()()222n n n k k k k k h h T f x f x f a f x f b --+==⎡⎤=+=++⎢⎥⎣⎦∑∑误差关于h 的函数：()()212n b a R fh f η-''=-页脚内容2复合辛普森公式：()()()()111/201426n n n k k k k h S f a f x f x f b --+==⎡⎤=+++⎢⎥⎣⎦∑∑误差关于h 的函数：()()441802n n b a h R f I S f η-⎛⎫=-=- ⎪⎝⎭1.2龙贝格求积算法：龙贝格求积公式是梯形法的递推化，也称为逐次分半加速法，它是在梯形公式、辛普森公式和柯特斯公式之间的关系的基础上，构造出一种计算积分的方法，同时它有在不断增加计算量的前提下提高误差的精度的特点。
计算过程如下：（1）取0,k h b a ==-，求：()()()[]()00.,.2hT f a f b k a b =+→⎡⎤⎣⎦令k 1记为区间的二分次数 （2）求梯形值02k b a T -⎛⎫⎪⎝⎭即按递推公式12102122n n n k k h T T f x -+=⎛⎫=+ ⎪⎝⎭∑计算0k T .（3）求加速值，按公式()()()111444141m m k k k mm m m m T T T +--=---逐个求出T 表的地k 行其余各元素()()1,2,,k j j T j k -=（4）若()()001k k T T ε--<（预先给定的精度），则终止计算，并取()()0;1k T I k k ≈+→否则令转（2）继续计算。
数值分析报告(第五版)计算实习的题目第三章
数值分析计算实习题第三章第二次作业：题一：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轴')。
数值分析计算实习题目一
数值分析计算实习题目一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 ====⋅=≤≤∏。
《数值分析》课程设计—16题
《数值分析》课程设计—作业实验一1.1水手、猴子和椰子问题：五个水手带了一只猴子来到南太平洋的一个荒岛上，发现那里有一大堆椰子。
由于旅途的颠簸，大家都很疲惫，很快就入睡了。
第一个水手醒来后，把椰子平分成五堆，将多余的一只给了猴子，他私藏了一堆后便又去睡了。
第二、第三、第四、第五个水手也陆续起来，和第一个水手一样，把椰子分成五堆，恰多一只猴子，私藏一堆，再去入睡，天亮以后，大家把余下的椰子重新等分成五堆，每人分一堆，正好余一只再给猴子，试问原先共有几只椰子？试分析椰子数目的变化规律，利用逆向递推的方法求解这一问题。
解：一、问题分析：对于本题，比较简单，我们只需要判断原来椰子的个数及每个人私藏了一份之后剩下的是否能被5除余1，直到最后分完。
对于第一个程序，n取2000；对于第二个程序，n取20001，就能得到我们想要的结果，即原先一共有15621个椰子，最终平均每人得4092个椰子。
1.2 当0,1,2,,100n =时，选择稳定的算法计算积分10d 10nxn xe I x e --=+⎰. 解：一、问题分析：由10d 10nxn xe I x e --=+⎰知： 1101001==+⎰dx I I 以及： )1(11010101010)1(1nnx x nx x n n n e ndx e dx e e e I I ----+-+-==++=+⎰⎰ 得递推关系： ⎪⎩⎪⎨⎧--=-=-+n nn I e n I I I 10)1(1101101， 但是通过仔细观察就能知道上述递推公式每一步都将误差放大十倍，即使初始误差很小，但是误差的传播会逐步扩大，也就是说用它构造的算法是不稳定的，因此我们改进上述递推公式（算法）如下：⎪⎪⎩⎪⎪⎨⎧--=-=+-))1(1(101)1(101110n n n I e n I I I通过比较不难得出该误差是逐步缩小的，即算法是稳定的。
二、问题求解：为了利用上面稳定的算法，需要我们估计初值100I 的值。
数值分析计算实习题
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]上三次样条插值更精确。
数值分析上机实验——数值积分
实验报告哈尔滨工程大学教务处制实验三 数值积分一．数值积分的基本思想1.复合梯形公式：Tn=++)()([2b f a f h2∑-=11)](n k xk f ；2.复合辛普森公式：Sn=6h[f(a)+f(b)+2∑-=11)](n k xk f +4∑-=+1)2/1(n k x f ]；以上两种算法都是将a-b 之间分成多个小区间（n ）,则h=(b-a)/n,x k =a+kh,x k+1/2=a+(k+1/2)h,利用梯形求积根据两公式便可。
3.龙贝格算法：在指定区间内将步长依次二分的过程中运用如下公式（1）Sn=34T2n-31Tn（2）Cn=1516S2n-151Sn（3）Rn=6364C2n-631Cn4T)(k m=144-m m T )1(1+-k m - 141-mT )(1k m -，k = 1，2，… 二.实验题目及实验目的(第4章计算实习题第1题)用不同数值方法计算积分xdx x ln 1⎰= -94。
（1）取不同的步长h 。
分别用复合梯形及复合辛普森求积计算积分，给出误差中关于h 的函数，并与积分精确值比较两个公式的精度，是否存在一个最小的h ，使得精度不能再被改善？（2）用龙贝格求积计算完成问题（1）。
（3）用自适应辛普森积分，使其精度达到104-。
三．实验手段：指操作环境和平台:win7系统下MATLAB R2009a程序语言：一种类似C 语言的程序语言，但比C 语言要宽松得多，非常方便。
四.程序①复合梯形求积程序function t=TiXing_quad(a,b,.h) format longx=a:h:b;y=sqrt(x).*log(x);y(1)=0;t=0;for k=1:(b-a)/h,t= t+y(k)+y(k+1);endt=t*h/2;②复合辛普森求积程序function s=Simpson_quad(a,b,h) format longx=a:h:b;y=sqrt(x).*log(x);z=sqrt(x+h/2).*log(x+h/2);y(1)=0;s=0;for k=1:(b-a)/h,s= s+y(k)+y(k+1)+4*z(k);ends=s*h./6;③龙贝格求积程序function [q,R]=Romberg(a,b,eps) h=b-a;R(1,1)=h*(0+sqrt(b).*log(b))/2; M=1;J=0;err=1;while err>epsJ=J+1;h=h/2;S=0;for p=1:Mx=a+h*(2*p-1);S=S+sqrt(x).*log(x);endR(J+1,1)=R(J,1)/2+h*S;M=2*M;for k=1:JR(J+1,k+1)=R(J+1,k)+(R(J+1,k)-R(J,k))/(4^k-1);enderr=abs(R(J+1,J)-R(J+1,J+1));endq=R(J+1,J+1);控制台输入代码：（1）>> a=0;>> b=1;>> h=0.1;>> t=TiXing_quad(a,b,h)>> s=Simpson_quad(a,b,h)>> h=0.01;>> t=TiXing_quad(a,b,h)>> s=Simpson_quad(a,b,h)>> h=0.001;>> t=TiXing_quad(a,b,h)>> s=Simpson_quad(a,b,h)（2）>> a=0;>> b=1;>> eps=10^-8;>> [quad,R]=Romberg(a,b,eps)（3）>> a=0;>> b=1;>> eps=10^-4;>> q=ZiShiYingSimpson('sqrt(x).*log(x)',a,b,eps) 五．实验结果比较与分析（1）h = 0.1时h = 0.01时h = 0.001时由结果（1）可知对于同一步长h，复合辛普森法求积分精度明显比复合梯形法求积的精度要高，且当步长取不同值时即h越小时，积分精度越高。
数值分析计算实习题列主元高斯消去法解线性方程组
数值分析计算实习题第5章解线性方程组的直接方法【选题列主元高斯消去法解线性方程组。
书上的计算实习题1、2、3都要求用列主元高斯消去法解线性方程组，所以考虑写一个普适的程序来实现。
对于线性方程组Ax二b,程序允许用户从文件读入矩阵数据或直接在屏幕输入数据。
文件输入格式要求:(1)第一行为一个整数n (2<=n<=100)，表示矩阵阶数。
(2)第2~n+l行为矩阵A各行列的值。
(3)第n+2~n+n+2行为矩阵b各行的值。
屏幕输入：按提示输入各个数据。
输出：A. b、det(A).列主元高斯消去计算过程、解向量X。
【算法说明】设有线性方程组Ax=b,其中设A为非奇异矩阵。
方程组的增广矩阵为«12«21[Nb] =第1步（k=l ）：首先在A的第一列中选取绝对值最大的元素®I,作为第一步的主元素:«|| H0然后交换（A, b）的第1行与第I行元素，再进行消元计算。
设列主元素消去法已经完成第1步到第k・l步的按列选主元，交换两行，消元计算得到与原方程组等价的方程组 A（k）x=b（k）4? …4;)…唸)•忒••輕■[A.b]T[A ⑹,b")] = ••■咲■■■■■* *■〃伏)・• - %■第k步计算如下: 对于 k=l, 2, •…，0-1（1）按列选主元：即确定t使（2）如果tHk,则交换［A, b］第t行与第k行元素。
（3）消元计算54* J 叫=一鱼(=^ + 1,…,H)% 吗 <-«y + 〃如伽 (fJ = R + l,…/)b- <-勺+加汝仇， (i = /c + l,…,《)消元乘数mik 满足:n (%-D 内)X1 < ------ -- ---- 9(j = « 一 1,«一2■…J)tk M 1,(,=斤 +1, •••,«)fet e(4)回代求解【程序】/*【普适列主元消去法解线性方程组】对于线性方程组：Ax=b 输入：［选择屏幕直接输入］1.A的行阶数n(l <= 11<= 100)2.A的值3.b的值［选择读取文件1文件名(和主程序同级文件夹下)输出：I.A2・b3・ det(A)4解向疑X */#inciude <stdio.h>#include <stdlibJi>#include <niath.h> double A[1051(J05LA_B[I05][105Lb[105].x[105];double det A:int n.mark = 1;〃读入数据void input(){int ij:char ch[20],name[ 100];HLE *f;printf("\iv-An是否从文件读取数据(Y/N):”); scanf(”％st&ch);if(ch[0] = Y II ch[0] = y)( prinif(“请输入文件名(包括扩展需)：")； scanf("%s".name):f = fop cn(namc「T');fscanfCf/'%d'\&n);ford = 0:i < n;i 卄) for(j = 0;j < nJ ++) fscanf(f,'*%ir\&A[i]|j))：for(i = 0:i < n;i卄)fscanf(「％F・&b[i]);else{prin氓”请输入A的阶数：”)；scanf( '%d %d\&n)： prinifC请输入A的值：”)；for(i = 0:i V n:i ++)for(j = 0;j < n:j ++) scanf("%lf\&A[i]U]);phnifC请输入b的值：”)；for(i = 0:i < n;i 卄)scanf(''%lf\&b[i])：〃讣算行列式的值double det{double s[105](105] jni m){int z.jkdouble b[105][105Klotal = 0.r; /*b[Nl[N]用于存放，在矩阵s[Nl[N冲元素s[0]的余子式灯if(m>2){for(z = 0：z V m:z++){for(j = 0;j V m ・ 1 j ++)for(k = 0:k vn卜l;k ++)if(k >= z)bUKk] = sU+l](k+l];elsebLi][k] = s[j+l][k];if(z % 2==0)r=s[0)⑵ * dcKb.m - 1)：/*递归调用制elser= (-1) * s[0](z] * det(b.m -1);total = total + r;else if(m == 2)total = s[0][0] *s(l][l]-s[0](l] *s[l][01：else if(m == 1)total =s[0][0];return total;// 输出A^llb和dcl(A)void ouipui_l(){ int i j;primlTAW);for(i = 0;i < n:i ++){ for(j = 0:j < n:j ++)p rintf('-%15.4f\A[i]lj]);prinif(W);prinlf(5b = \rf);for{i = 0;i < n:i ++)prinif("%15.4l\iV\b[i]):printf("\ndet(A) = %・4f\n”・dclA);//主il•算函数void couni_x(){int ij,k:int max; double tmpjnik;〃构造增广矩阵for{i = 0;i<n:i++){for。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
数值分析——数值积分实习题
管理科学与工程学院 学号：1120140500 姓名：彭洋洋 一、计算实习题
1.
用不同数值方法计算积分：04
9
xdx =-⎰.
（1）取不同的步长h ，分别用复合梯形及复合辛普森求积计算积分，给出误差中关于h 的函数，并与积分精确值比较两个公式的精度，是否存在一个最小的h ，使得精度不能再被改善？ （2）用龙贝格求积计算完成问题（1） （3）用自适应辛普森积分，使其精度达到10-4
解答：（1）取不同的步长，采用不同的公式，比较精度过程如下： 1.1 复合梯形公式及复合辛普森公式求解
复合梯形公式：1
1
*[()2()()]2n n k k h
T f a f x f b -==++∑
误差关于h 的函数：2
(2)()**()12
n a b R f h f ξ-=
复合辛普森公式：11
1/201
*[()4()2()()]6n n n k k k k h
S f a f x f x f b --+===+++∑∑
误差关于h 的函数：4(4)()*(/2)*()180
n a b
R f h f η-=
1.2 复合梯形公式及复合辛普森公式Matlab 程序
（2）用龙贝格求积计算完成问题（1） 2.1 龙贝格求积算法
龙贝格求积公式也称为逐次分半加速法。
它是在梯形公式、辛普森公式和柯特斯公式之间的关系的基础上，构造出一种加速计算积分的方法。
作为一种外推算法，它在不增加计算量的前提下提高了误差的精度。
24133n n n S T T =- 21611515n n n C S S =- 26416363
n n n R C C =-
1221/20
1()22n n n k k h T T f x -+==+∑ ()(1)()11
(4*)/(41)k m k k m
m m m T T T +--=-- 1,2,...k = 2.2 龙贝格求积Matlab 程序
画图程序设计
①得到关于n各种公式求积的图表如下：
对于梯形公式、辛普森公式n代表份数，龙贝格公式n表示从1开始的序列号
②关于步长h 的各种公式求积的图表如下
其中龙贝格序列步长()/2k h b a =-：
观察两幅图表h 越小，精度越高。
当h 足够小时，精度没有必要再被改善。
由结果（1）、（2）可知复合辛普森法求积分精度明显比复合梯形法求积的精度要高，且当步长取不同值时即n 越大、h 越小时，积分精度越高。
实验结果说明不存在一个最小的h ，使得精度不能再被改善。
又两个相应的关于h 的误差（余项）2(2)()**()12n a b R f h f ξ-=
和4(4)()*(/2)*()180
n a b
R f h f η-= 其中,(a,b)ξη∈，可知h 愈小，余项愈小，从而积分精度越高。
附（部分表格）：
（3）用自适应辛普森积分，使其精度达到10-4
自适应积分法是一种计较经济而且快速的求积分的方法，它能自动地在被积函数变化剧烈的地方增多节点，而在被积函数变化平缓的地方减少节点，是一种不均匀区间的积分方法。
求解结果如下
>> F=inline('sqrt(x).*log(x)');
>> s=quad(F,0,1,1e-4) s = -0.4440。