数值分析上机实验6
例1．世界人口数据拟合问题：据统计，六十年代世界人口数据如下(单位：亿)
根据表中数据，预测公元2000年时的世界人口。
问题分析与数学模型
设人口总数为N(t)，根据人口理论的马尔萨斯模型，采用指数函数
N(t) = e a + b t
=+，令
对数据进行拟合。为了计算方便，将上式两边同取对数，得ln N a bt
y = ln N或N = e y
变换后的拟合函数为
y(t) = a + b t
根据表中数据及等式 k k( 1，2，……，9)可列出关于两个未知数、b的9个方程的超定方程组（方程数多于未知数个数的方程组）
a + t j b = y j(j= 1，2， (9)
可用最小二乘法求解。
算法与数学模型求解
算法如下：
第一步：输入人口数据，并计算所有人口数据的对数值；
第二步：建立超定方程组的系数矩阵，并计算对应的正规方程组的系数矩阵和右端向量；第三步：求解超定方程组并输出结果：a，b；
第四步：利用数据结果构造指数函数计算2000年人口近似值N(2000)，结束。
MATLAB程序
t=1960:1968;t0=2000;
N=[29.72 30.61 31.51 32.13 32.34 32.85 33.56 34.20 34.83];
y=log(N);
A=[ ones(9,1), t' ];
d=A\ y' ;a=d(1),b=d(2)
N0=exp(a+b*t0)
x=1960:2001;yy=exp(a+b*x);
plot(x,yy,t,N,'o',2000,N0,'o')
计算结果为
a =－33.0383，
b = 0.0186
N (2000) = 63.2336
所以取五位有效数，可得人口数据的指数拟合函数
t e t N 0186.00383.33)(+-=
经计算得2000年人口预测值为：63.2336 (亿)。
例2．温度数据的三角函数拟合问题 洛杉矶郊区在11月8日的温度记录如下
在不长的时期内，气温的变化常以24小时为周期，考虑用Fourier 级数的部分和（有限项）做拟合函数。即，求最小二乘曲线：
∑=++=n
k k k x k b x k a a x 1
0)]sin()cos([)(ωω?
其中，ω = 2π / 24。例如，当n=2时拟合函数为
?(x )= a 0 + a 1 cos(ωx ) + b 1 sin(ωx ) + a 2 cos(2ωx ) + b 2 sin(2ωx )
对不同的n ，确定拟合函数中的各系数。绘出最小二乘曲线与离散数据点，并计算出拟合函
数的残差2-范数。 算法分析：
以n=2时的拟合函数为对象作算法分析。24小时的温度记录可列为数表如下
将24个数据点代入拟合函数得超定方程组
??
?????
?
????????=???????????????????????????
??
???2421221102424
2424
22
2
2111
12sin 2cos sin cos 12sin 2cos sin cos 12sin 2cos sin cos 1y y y b a b a a x x x x x x x x
x x x x
ωωωωωωωωωωωω 可以证明方程组的系数矩阵列向量组是正交向量组，于是由最二乘法所推出的正规方程组
系数矩阵是对角矩阵。所以原方程组的最小二乘解为
∑==24
1
0241k k y a
∑∑===24
12
24
11cos /cos k k k k k x y x a ωω，∑∑===24
1224
122cos /2cos k k k k k x y x a ωω
∑∑===24
1
224
1
1sin /sin k k k k k x y x b ωω， ∑∑===24
1
224
1
22sin /2sin k k k k k x y x b ωω
MATLAB 程序（运行程序时需输入参数n ）： n=input('input n=: '); w=2*pi/24;x=[1:24]';
y=[66;66;65;64;63;63;62;61;60;60;59;58; 58;58;58;58;57;57;57;58;60;64;67;68]; a0=sum(y)/24; for k=1:n
ck=cos(k*w*x);sk=sin(k*w*x);
a(k)=(ck'*y)/(ck'*ck); b(k)=(sk'*y)/(sk'*sk); end yy=a0; for k=1:n
yy=yy+a(k)*cos(k*w*x)+b(k)*sin(k*w*x); end
plot(x,y,'x',x,yy) r=norm(yy-y)
3．切比雪夫多项式的前两项为：T 0(x ) = 1，T 1(x ) = x ，对于n ≥2，有递推公式
T n+1(x ) = 2xT n (x ) – T n – 1(x )
当x ∈[ – 1，1 ] 时，利用递推公式，计算并绘出 T 0(x )，T 1(x )，T 2(x )，T 3(x )，T 4(x )的
函数图形
MATLAB 程序如下：
x=-1:.05:1;
T0=ones(size(x)); T1=x;
plot(x,T0,'b',x,T1,'b'); hold on for k=2:4
T=2*x.*T1-T0; plot(x,T)
T0=T1;T1=T; end axis off
4．1912年，伯恩斯坦给出了关于多项式一致逼近连续函数的构造性证明，提出了著名的伯恩斯坦多项式，设 f (x )在区间 [0，1]上连续，他的多项式为
B x f k
n C x x n n k n k k k n
()()()=--=∑10
试利用组合数的递推公式 1
11---+=k n k n k n C C C ，设计一个计算n 次伯恩斯坦多项式函数值的
算法。并对函数 f (x ) = sin x 给以验证。 MATLAB 程序如下 n=input('input n='); x=[0:n]/n; f=sin(x*pi); for i=1:n+1 y=f;t=x(i); for k=n:-1:1 for j=1:k
y(j)=t*y(j)+(1-t)*y(j+1); end end
p(i)=y(1);
end
max(abs(f-p))
plot(x,f,'b',x,p,'o',x,p,'r')
运行四次程序，分别输入n=10，n=20，n=30，n=40得下面图形
5．Bezier曲线是法国雷诺汽车公司的工程师Bezier于1971年提出了一种新的参数曲线表示法。这种方法可以交互式地确定一组控制多边形顶点以获得所需要的曲线形式。设曲线参数方程
x = x(t)，y = y(t)
如果给定控制多边形顶点P0，P1，…，P m的坐标
(x0，y0)，(x1，y1)，……，(x m，y m)
则相应的Bezier多项式由下式定义
∑=
--
=
m
k
k
k
m
k
k
m
x
t
t
C
t
x
) 1(
)(，∑
=
--
=
m
k
k
k
m
k
k
m
y
t
t
C
t
y
) 1(
)(
用矢量函数的形式来表示平面曲线则有数学表达式
∑=--=m
k k k m k
k m P t t C t P 0)1()(
三次Bezier 曲线数学表达式为
P t t P t t P t t P t P t ()()()()=-+-+-+≤≤1313101302122
3
3
试给定四个平面点绘制飞机机翼剖面轮廓曲线 function z=bez3(p) t=[0:.05:1]';t1=1-t;
z=[t1.^3 3*t.*t1.^2 3*t1.*t.^2 t.^3]*p; plot(p(:,1),p(:,2),z(:,1),z(:,2))
首先录入这一函数文件，文件名为bez3.m 。为了调用这一函数，必须先定义多边折线的结点数组，在MA TLAB 环境下键入： p=[0 0;0 0.4;0.5 1;2 0] bez3(p)
计算机将绘制出对应的曲线。
一、数值试验
1. SARS 的传播及预防问题
非典的爆发和蔓延给我国的经济发展和人民生活带来了很大的影响，下表给出了北京市当年4月份到6月份的疫情数据,通过拟合确诊的累积病人曲线,若延后5天采取严格的预防措施,对疫情的传播所生成的影响做出估计. 日期
已确诊病例累积
现有疑似病例
死亡累积
治愈出院累积 4月20日 297 402 18 33 4月30日 1584 1408 75 90 5月1日 1640 1415 82 100 5月10日 1988 1397 116 175 5月20日 2189 1225 150 395 5月30日
2309
706
176
1006
6月1日 2319 739 181 1124 6月10日
2394
351
184 1747 6月20日 2439 3
191
2189
(1) 根据已有数据，用MATLAB 绘出确诊累积病人离散图
(2) 取拟合曲线的拟合函数为如下非线性函数
x
b a y +=1 试确定拟合函数中的参数：a ，b ，并推测五天后累积病人数量。
2．蠓虫分类问题
生物学家试图对两类蠓虫(Af 与 Apf)进行鉴别,依据的资料是蠓虫的触角和翅膀的长度,已经测得9只Af 和6只Apf 的数据(触角长度用x 表示,翅膀长度用y 表示)
现需要解决三个问题：
（1）如何凭借原始资料(15 对数据,被称之为学习样本)制定一种方法区分两类蠓虫； （2）依据确立的方法,对题目提供的三个样本：(1.24,1.80)，(1.28,1.84)，(1.40,2.04)加以识别；
（3）设Af 是宝贵的传粉益虫，Apf 是某种疾病的载体，是否应该修改分类方法。
问题分析：首先画出15对数据的散点图，其中，Af 用 * 标记，Apf 用 × 标记。观察图形，可以发现，Af 的点集中在图中右下角，而 Apf 的点集中在图中左上角。应该存在一条直线 L 位于两类点之间， 作为 Af 和 Apf 分界线，这条直线 L 的确定应依据问题所给的数据，即学习样本。设这条直线的方程为
w 1 x + w 2 y + w 0 = 0
对于平面上任意一点 P (x ，y )，如果该点在直线上，将其坐标代入直线方程则使方程成为恒等式，即使方程左端恒为零；如果点P (x ，y )不在直线上，将其坐标代入直线方程，则方程左端不为零。由于 Af 和 Apf 的散点都不在所求的直线上，故将问题所提供的数据代入直线方程左端应该得到表达式的值大于零或者小于零两种不同的结果。
为了建立判别准则，利用直线方程左端表达式定义判别函数g (x ，y )，使得当平面上某点P (x ，y )属于 Af 类时，g (x ，y )>0，否则 g (x ，y )<0。
引入学习机制，在学习过程中将两种不同的状态，以“1”和“－1”表示。当P (x ，y )属于 Af 类时，g (x ，y )=1，否则 g (x ，y )=－1。取
g (x ，y )= w 1 x + w 2 y + w 0
其中，w 1，w 2，w 0为三个待定系数，由所给数据形成约束条件，即线性方程组：
w x w y w j w x w y w j j j j j 120120112911015++==++=-=??
??,
(,,,),
(,,)
这是包括三个未知数共15个方程的超定方程组，可以求方程组的最小二乘解。
实验要求：
（1）根据上面分析写出对应的正规方程组并求解；
（2）确定分类边界直线的方程。由所给数据用判别函数判别三个新蠓虫的类属； （3）如果Apf 是某种疾病的载体，应该修改数学模型确定分类方法。 3．海浪高度数据的拟合问题
海洋水文观测站记录了某海域每天24小时海浪潮高度数据（相对于海堤上的零标尺记号）。工作人员每小时做一次记录，下面是从某年12月1日早上零点开始到12月2日晚上23点结束两天内完整的数据记录。能否根据这些数据预测三天以后（即12月5日）的海浪高度数据。
二、方法应用与编程练习
1．龙格（Range ）反例是插值方法中一个著名的反例。比较多项式插值和多项式拟合的方法的异同，借用MA TLAB 中数据拟合命令做龙格反例的数值实验。 2根据马尔萨斯人口模型，设利用数据拟合方法确定指数函数，并预测2000年到2008年九年的我国人口数量。
3．勒让德多项式的前两项为：p 0(x ) = 1，p 1(x ) = x ，对于n ≥2，有递推公式
p n+1(x ) = [(2n+1) x p n (
x ) – n p n – 1(x )]/(n+1)
对x ∈[ – 1，1 ] ，绘制勒让德多项式p 0(x )，p 1(x )，p 2(x )，p 3(x )，p 4(x )图形
4
由
?
cos 1e p
r -=
其中，p 为参数，e 为偏心率，由给定的观察值用最小二乘方法拟合轨道曲线，求出p 和e ，并给出平方误差。
令k=1/r，处理数据。解超定方程组确定p和e
r=[2.70 2.00 1.61 1.20 1.02];
fai=[48 67 83 108 126]*pi/180;
n=length(r);b=ones(n,1);
G=[1./r' cos(fai')];
X=G\b;
p=X(1);e=X(2);
fa=pi/4:.1:3.5*pi/2;
ro=p./(1-e*cos(fa));
polar(fai,r,'ro'),hold on
polar(fa,ro)
5．斐波拉奇数列中整数的高精度计算问题。对n=100和n=200，分别计算斐波拉奇数列的第n项和前n项和（要求全部用整数表示）。
6．多边形面积的计算问题。任意一个n边形都可以分割成为n-2个三角形，分别计算出这些三角形的面积，将所有三角形面积和作为n边形面积是一种解决问题的有效方法。
输入数据的第一行有一个正整数n，表示多边形顶点数目（0≤n≤1000）。接下来有n行数据，每一行有一对实数x i，y i 用来表示多边形顶点的平面坐标（假定这些顶点确实构成多边形）。对输入数据计算出多边形面积S。
Sample Imput 5
0 0
0 1
0.5
0.5
1 1
1 0
Output for the Sample Input 0.75
数值分析上机题 第一章 17.（上机题）舍入误差与有效数 设∑=-= N j N j S 2 2 11 ，其精确值为）111-23（21+-N N 。 （1）编制按从大到小的顺序1 -1 ···1-311-21222N S N +++=，计算N S 的通用 程序； （2）编制按从小到大的顺序1 21 ···1)1(111 222-++--+ -=N N S N ，计算N S 的通用程序； （3）按两种顺序分别计算210S ,410S ,610S ,并指出有效位数（编制程序时用单精度）； （4）通过本上机题，你明白了什么？ 解： 程序： （1）从大到小的顺序计算1 -1 ···1-311-21222N S N +++= ： function sn1=fromlarge(n) %从大到小计算sn1 format long ; sn1=single(0); for m=2:1:n sn1=sn1+1/(m^2-1); end end （2）从小到大计算1 21 ···1)1(111 2 22 -++--+-= N N S N function sn2=fromsmall(n) %从小到大计算sn2 format long ; sn2=single(0); for m=n:-1:2 sn2=sn2+1/(m^2-1); end end （3） 总的编程程序为： function p203()
clear all format long; n=input('please enter a number as the n:') sn=1/2*(3/2-1/n-1/(n+1));%精确值为sn fprintf('精确值为%f\n',sn); sn1=fromlarge(n); fprintf('从大到小计算的值为%f\n',sn1); sn2=fromsmall(n); fprintf('从小到大计算的值为%f\n',sn2); function sn1=fromlarge(n) %从大到小计算sn1 format long; sn1=single(0); for m=2:1:n sn1=sn1+1/(m^2-1); end end function sn2=fromsmall(n) %从小到大计算sn2 format long; sn2=single(0); for m=n:-1:2 sn2=sn2+1/(m^2-1); end end end 运行结果：
实验一 误差分析 实验（病态问题） 实验目的：算法有“优”与“劣”之分，问题也有“好”与“坏”之别。对数值方法的研究而言，所谓坏问题就是问题本身对扰动敏感者，反之属于好问题。通过本实验可获得一个初步体会。 数值分析的大部分研究课题中，如线性代数方程组、矩阵特征值问题、非线性方程及方程组等都存在病态的问题。病态问题要通过研究和构造特殊的算法来解决，当然一般要付出一些代价（如耗用更多的机器时间、占用更多的存储空间等）。 问题提出：考虑一个高次的代数多项式 )1.1() ()20()2)(1()(20 1∏=-=---=k k x x x x x p 显然该多项式的全部根为1,2,…,20共计20个，且每个根都是单重的。现考虑该多项式的一个扰动 )2.1(0 )(19=+x x p ε 其中ε是一个非常小的数。这相当于是对（）中19x 的系数作一个小的扰动。我们希望比较（）和（）根的差别，从而分析方程（）的解对扰动的敏感性。 实验内容：为了实现方便，我们先介绍两个Matlab 函数：“roots ”和“poly ”。 roots(a)u = 其中若变量a 存储n+1维的向量，则该函数的输出u 为一个n 维的向量。设a 的元素依次为121,,,+n a a a ，则输出u 的各分量是多项式方程 01121=+++++-n n n n a x a x a x a 的全部根；而函数 poly(v)b =
的输出b 是一个n+1维变量，它是以n 维变量v 的各分量为根的多项式的系数。可见“roots ”和“poly ”是两个互逆的运算函数。 ;000000001.0=ess );21,1(zeros ve = ;)2(ess ve = ))20:1((ve poly roots + 上述简单的Matlab 程序便得到（）的全部根，程序中的“ess ”即是（）中的ε。 实验要求： （1）选择充分小的ess ，反复进行上述实验，记录结果的变化并分析它们。 如果扰动项的系数ε很小，我们自然感觉（）和（）的解应当相差很小。计算中你有什么出乎意料的发现表明有些解关于如此的扰动敏感性如何 （2）将方程（）中的扰动项改成18x ε或其它形式，实验中又有怎样的现象 出现 （3）（选作部分）请从理论上分析产生这一问题的根源。注意我们可以将 方程（）写成展开的形式， ) 3.1(0 ),(1920=+-= x x x p αα 同时将方程的解x 看成是系数α的函数，考察方程的某个解关于α的扰动是否敏感，与研究它关于α的导数的大小有何关系为什么你发现了什么现象，哪些根关于α的变化更敏感 思考题一：（上述实验的改进） 在上述实验中我们会发现用roots 函数求解多项式方程的精度不高，为此你可以考虑用符号函数solve 来提高解的精确度，这需要用到将多项式转换为符号多项式的函数poly2sym,函数的具体使用方法可参考Matlab 的帮助。
实验一误差分析 实验1.1（病态问题） 实验目的：算法有“优”与“劣”之分，问题也有“好”与“坏”之别。对数值方法的研究而言，所谓坏问题就是问题本身对扰动敏感者，反之属于好问题。通过本实验可获得一个初步体会。 数值分析的大部分研究课题中，如线性代数方程组、矩阵特征值问题、非线性方程及方程组等都存在病态的问题。病态问题要通过研究和构造特殊的算法来解决，当然一般要付出一些代价（如耗用更多的机器时间、占用更多的存储空间等）。 问题提出：考虑一个高次的代数多项式 显然该多项式的全部根为1,2,…,20共计20个，且每个根都是单重的。现考虑该多项式的一个扰动 其中ε（1.1）和（1.221,,,a a 的输出b ”和“poly ε。 （1（2 （3）写成展 关于α solve 来提高解的精确度，这需要用到将多项式转换为符号多项式的函数poly2sym,函数的具体使用方法可参考Matlab 的帮助。 实验过程： 程序： a=poly(1:20); rr=roots(a); forn=2:21 n form=1:9 ess=10^(-6-m);
ve=zeros(1,21); ve(n)=ess; r=roots(a+ve); -6-m s=max(abs(r-rr)) end end 利用符号函数：（思考题一）a=poly(1:20); y=poly2sym(a); rr=solve(y) n
很容易的得出对一个多次的代数多项式的其中某一项进行很小的扰动，对其多项式的根会有一定的扰动的，所以对于这类病态问题可以借助于MATLAB来进行问题的分析。 学号：06450210 姓名：万轩 实验二插值法
数值分析实习报告 姓名：gestepoA 学号：201******* 班级：***班
序言 随着计算机技术的迅速发展，数值分析在工程技术领域中的应用越来越广泛，并且成为数学与计算机之间的桥梁。要解决工程问题，往往需要处理很多数学模型，不仅要研究各种数学问题的数值解法，同时也要分析所用的数值解法在理论上的合理性，如解法所产生的误差能否满足精度要求：解法是否稳定、是否收敛及熟练的速度等。而且还能减少大量的人工计算。 由于工程实际中所遇到的数学模型求解过程迭代次数很多，计算量很大，所以需要借助如MATLAB，C++，VB，JAVA的辅助软件来解决，得到一个满足误差限的解。本文所计算题目，均采用MATLAB进行编程，MATLAB被称为第四代计算机语言，利用其丰富的函数资源，使编程人员从繁琐的程序代码中解放出来MATLAB最突出的特点就是简洁，它用更直观的、符合人们思维习惯的代码。它具有以下优点： 1友好的工作平台和编程环境。MATLAB界面精致，人机交互性强，操作简单。 2简单易用的程序语言。MATLAB是一个高级的矩阵/阵列语言，包含控制语言、函数、数据结构，具有输入、输出和面向对象编程特点。用户可以在命令窗口中将输入语句与执行命令同步，也可以先编好一个较大的复杂的应用程序（M 文件）后再一起运行。 3强大的科学计算机数据处理能力。包含大量计算算法的集合，拥有600多个工程中要用到的数学运算函数。 4出色的图像处理功能，可以方便地输出二维图像，便于我们绘制函数图像。
目录 1 第一题 (4) 1.1 实验目的 (4) 1.2 实验原理和方法 (4) 1.3 实验结果 (5) 1.3.1 最佳平方逼近法 (5) 1.3.2 拉格朗日插值法 (7) 1.3.3 对比 (8) 2 第二题 (9) 2.1实验目的 (9) 2.2 实验原理和方法 (10) 2.3 实验结果 (10) 2.3.1 第一问 (10) 2.3.2 第二问 (11) 2.3.3 第三问 (11) 3 第三题 (12) 3.1实验目的 (12) 3.2 实验原理和方法 (12) 3.3 实验结果 (12) 4 MATLAB程序 (14)
数值分析实验指导 2012年8月
实验一 误差分析 实验1.1（病态问题） 实验目的：算法有“优”与“劣”之分，问题也有“好”与“坏”之别。对数值方法的研究而言，所谓坏问题就是问题本身对扰动敏感者，反之属于好问题。通过本实验可获得一个初步体会。 数值分析的大部分研究课题中，如线性代数方程组、矩阵特征值问题、非线性方程及方程组等都存在病态的问题。病态问题要通过研究和构造特殊的算法来解决，当然一般要付出一些代价（如耗用更多的机器时间、占用更多的存储空间等）。 问题提出：考虑一个高次的代数多项式 )1.1() ()20()2)(1()(20 1∏=-=---=k k x x x x x p 显然该多项式的全部根为1,2,…,20共计20个，且每个根都是单重的。现考虑该多项式的一个扰动 )2.1(0 )(19=+x x p ε 其中ε是一个非常小的数。这相当于是对（1.1）中19x 的系数作一个小的扰动。我们希望比较（1.1）和（1.2）根的差别，从而分析方程（1.1）的解对扰动的敏感性。 实验内容：为了实现方便，我们先介绍两个MATLAB 函数：“roots ”和“poly ”。 roots(a)u = 其中若变量a 存储n+1维的向量，则该函数的输出u 为一个n 维的向量。设a 的元素依次为121,,,+n a a a ，则输出u 的各分量是多项式方程 01121=+++++-n n n n a x a x a x a 的全部根；而函数 p o l y (v b = 的输出b 是一个n+1维向量，它是以n 维向量v 的各分量为根的多项式的系数。可见“roots ”和“poly ”是两个互逆的运算函数。 ) )20:1((; )2();21,1(;000000001.0ve poly roots ess ve zeros ve ess +=== 上述简单的MATLAB 程序便得到（1.2）的全部根，程序中的“ess ”即是（1.2）中的ε。
数值分析上机实验报告
《数值分析》上机实验报告 1.用Newton 法求方程 X 7-X 4+14=0 在（0.1，1.9）中的近似根（初始近似值取为区间端点，迭代6次或误差小于0.00001）。 1.1 理论依据： 设函数在有限区间[a ，b]上二阶导数存在，且满足条件 {}α?上的惟一解在区间平方收敛于方程所生的迭代序列 迭代过程由则对任意初始近似值达到的一个中使是其中上不变号 在区间],[0)(3,2,1,0,) (') ()(],,[x |))(),((|,|,)(||)(|.4;0)(.3],[)(.20 )()(.110......b a x f x k x f x f x x x Newton b a b f a f mir b a c x f a b c f x f b a x f b f x f k k k k k k ==- ==∈≤-≠>+ 令 )9.1()9.1(0)8(4233642)(0)16(71127)(0)9.1(,0)1.0(,1428)(3 2 2 5 333647>?''<-=-=''<-=-='<>+-=f f x x x x x f x x x x x f f f x x x f 故以1.9为起点 ?? ?? ? ='- =+9.1)()(01x x f x f x x k k k k 如此一次一次的迭代，逼近x 的真实根。当前后两个的差<=ε时，就认为求出了近似的根。本程序用Newton 法求代数方程(最高次数不大于10)在（a,b ）区间的根。
1.2 C语言程序原代码： #include
数值分析第一次实验报告 姓名： 学号： 实验1： 1. 实验项目的性质和任务 通过上机实验，使学生对病态问题、线性方程组求解和函数的数值逼近方法有一个初步理解。 2．教学内容和要求 1）对高阶多多项式 20 1()(1)(2)(20)()k p x x x x x k ==---=-∏ 编程求下面方程的解 19()0p x x ε+= 并绘图演示方程的解与扰动量ε的关系。（实验） 2）对2~20n =，生成对应的Hilbert 矩阵，计算矩阵的条件数；通过先确定解获得常向量b 的方法,确定方程组 n H x b = 最后，用矩阵分解方法求解方程组，并分析计算结果。（第三章，实验题4） 3）对函数 2 1()[1,1]125f x x x =∈-+ 的Chebyshev 点 (21)cos( ) 1,2,...,12(1) k k x k n n π -==++ 编程进行Lagrange 插值，并分析插值结果。（第四章 实验1）
项目涉及核心知识点 病态方程求解、矩阵分解和方程组求解、Lagrange插值。 重点与难点 算法设计和matlab编程。 1）a．实验方案： 先创建一个20*50的零矩阵X，然后利用Matlab中的roots（）和poly（）函数将50个不同的ess扰动值所产生的50个解向量分别存入X矩阵中。然后再将ess向量分别和X的20个行向量绘图。即可直观的看出充分小的扰动值会产生非常大的偏差。即证明了这个问题的病态性。 b．编写程序： >> X=zeros(20,50); >> ve=zeros(1,21); >> ess=linspace(0,,50);k=1; >> while k<=50 ve(2)=ess(k); X(1:20,k)=roots(poly(1:20)+ve); k=k+1; end >> m=1; >> while m<=20 figure(m),plot(ess,X(m,:));
数值分析上机实验报告 课程名称：数值分析上机实验 学院：机械工程学院专业：机械制造 姓名：张法光学号：2012021691 年级：12级任课教师：代新敏老师 2012年12月30日
一．已知A 与b 12.38412 2.115237 -1.061074 1.112336 -0.1135840.718719 1.742382 3.067813 -2.031743 2.11523719.141823 -3.125432 -1.012345 2.189736 1.563849 -0.784165 1.112348 3.123124 -1.061074 -3.125A =43215.567914 3.123848 2.031454 1.836742-1.056781 0.336993 -1.010103 1.112336 -1.012345 3.12384827.108437 4.101011-3.741856 2.101023 -0.71828 -0.037585 -0.113584 2.189736 2.031454 4.10101119.8979180.431637- 3.111223 2.121314 1.784317 0.718719 1.563849 1.836742 -3.741856 0.4316379.789365-0.103458 -1.103456 0.238417 1.742382 -0.784165 -1.056781 2.101023-3.111223-0.1034581 4.7138465 3.123789 -2.213474 3.067813 1.112348 0.336993-0.71828 2.121314-1.103456 3.12378930.719334 4.446782 -2.031743 3.123124 -1.010103-0.037585 1.7843170.238417-2.213474 4.44678240.00001[ 2.1874369 33.992318 -2 5.173417 0.84671695 1.784317 -8 6.612343 1.1101230 4.719345 -5.6784392]T B ????? ? ?? ? ? ???? ? ? ???? ? ? ????? ?=（2）用超松弛法求解Bx=b （取松弛因子ω=1.4，x （0）=0，迭代9次）。 （3）用列主元素消去法求解 Bx=b 。 解：（3）、用列主元素消去法求解Bx=b （一）、理论依据： 其基本思想是选取绝对值尽量大的元素作为主元素，进行行与列的交换，再进行回代，求出方程的解。 将方阵A 和向量b 写成C=（A b ）。将C 的第1列中第1行的元素与其下面的此列的元素逐一进行比较，找到最大的元素1j c ，将第j 行的元素与第1行的元素进行交换，然后通过行变换，将第1列中第2到第n 个元素都消成0。将变换后的矩阵(1)C 的第二列中第二行的元 素与其下面的此列的元素逐一进行比较，找到最大的元素(1) 2k c ，将第k 行的元素与第2行的 元素进行交换，然后通过行变换，将第2列中第3到第n 个元素都消成0。以此方法将矩阵的左下部分全都消成0。 （二）、计算程序： #include "math.h" #include "stdio.h" void main() { double u[9],x1[9],y[9],q[9],b1[9][10],x[9],a[9][9]={ {12.38412,2.115237,-1.061074,1.112336,-0.113584,0.718719,1.742382,3.067813,-2.031743 },
1.给定n 阶方程组Ax b =，其中 6186186186A ?? ? ? ?= ? ? ??? ，7151514b ?? ? ? ?= ? ? ??? 则方程组有解(1,1,,1)T x = 。对10n =和84n =，分别用Gauss 消去法和列主元消去法解方程组，并比较计算结果。 Gauss 消去法： Matlab 编程（建立GS.m 文件）： function x=GS(n) A=[];b=[]; for i=1:n-1 A(i,i)=6; A(i,i+1)=1; A(i+1,i)=8; b(i)=15; end A(n,n)=6;b(1)=7;b(n)=14;b=b'; for k=1:n-1 for i=k+1:n m(i,k)=A(i,k)/A(k,k); A(i,k:n)=A(i,k:n)-m(i,k)*A(k,k:n); b(i)=b(i)-m(i,k)*b(k); end end b(n)=b(n)/A(n,n); for i=n-1:-1:1 b(i)=(b(i)-sum(A(i,i+1:n).*b(i+1:n)'))/A(i,i); end clear x; x=b; disp( 'AX=b 的解x 是') end 计算结果： 在matlab 命令框里输出GS （10）得： >> GS(10) AX=b 的解x 是 ans = 1.0000 1.0000
1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 在matlab命令框里输出GS（84）得：>> GS(84) AX=b的解x是 ans = 1.0e+008 * 0.0000 … … … 0.0000 -0.0000 0.0000 -0.0000 0.0000 -0.0000 0.0000 -0.0000 0.0000 -0.0000 0.0000 -0.0000 0.0000 -0.0001 0.0002 -0.0003 0.0007 -0.0013 0.0026 -0.0052 0.0105 -0.0209 0.0419 -0.0836 0.1665 -0.3303
内容包括： 实验题目1：算法的数值稳定性实验 实验题目2：LU分解实验 实验题目3：三次样条插值外推样条实验 实验题目4：第二类Fredholm 积分方程实验实验题目5：M级显式R_K法
实验题目：算法的数值稳定性实验 实验内容：计算积分()1 0()d 1515n x I n x a x ==+? (n=1,2,…,20) 易得到下面递推公式 ()()1 1I n aI n n =--+ 并有估计式 ()() ()() 1 1 111I n a n a n << +++ 计算方法： 算法一:采用下面递推公式计算: ()()1 1I n aI n n =--+ ()1,2,,20 n = 取初值()116 0ln ln 15a I a +== 算法二: 采用下面递推公式计算: ()()111I n I n a n ??-= -+???? ()20,19,,1 n =
结果分析：（分析哪个好哪个不好，原因是什么） 我觉得算法二比较好, 原因一:根据式 ()() ()() 1 1 111I 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 εε-= ,…, ^ ^ 01 n n a εε=,尽管有初始误差^ 20ε,但随着计算的进程,这个误差的影响不断减小。 附：源程序：（把源程序附上） 算法一程序: >> format long >> a=15;I=log(16/15); for n=1:20 n I=-a*I+1/n end 算法二程序: >> format long >> a=15;I=31/10080; >> for n=20:-1:1 n I I=1/a*(-I+1/n); End
课题一：拉格朗日插值法 1.实验目的 1．学习和掌握拉格朗日插值多项式。 2.运用拉格朗日插值多项式进行计算。 2.实验过程 作出插值点（1.00，0.00），（-1.00，-3.00），（2.00，4.00）二、算法步骤 已知：某些点的坐标以及点数。 输入：条件点数以及这些点的坐标。 输出：根据给定的点求出其对应的拉格朗日插值多项式的值。 3.程序流程： （1）输入已知点的个数； （2）分别输入已知点的X坐标； （3）分别输入已知点的Y坐标； 程序如下： #include
插值算法*/ { int i,j; float *a,yy=0.0; /*a作为临时变量，记录拉格朗日插值多项*/ a=(float*)malloc(n*sizeof(float)); for(i=0;i<=n-1;i++) { a[i]=y[i]; for(j=0;j<=n-1;j++) if(j!=i) a[i]*=(xx-x[j])/(x[i]-x[j]); yy+=a[i]; } free(a); return yy; } int main() { int i; int n; float x[20],y[20],xx,yy; printf("Input n:");
scanf("%d",&n); if(n<=0) { printf("Error! The value of n must in (0,20)."); getch();return 1; } for(i=0;i<=n-1;i++) { printf("x[%d]:",i); scanf("%f",&x[i]); } printf("\n"); for(i=0;i<=n-1;i++) { printf("y[%d]:",i);scanf("%f",&y[i]); } printf("\n"); printf("Input xx:"); scanf("%f",&xx); yy=lagrange(x,y,xx,n); printf("x=%f,y=%f\n",xx,yy); getch(); } 举例如下：已知当x=1，-1,2时f(x)=0，-3,4，求f(1.5)的值。
实验报告一 题目：非线性方程求解 摘要：非线性方程的解析解通常很难给出，因此线性方程的数值解法就尤为重要。本实验采用两种常见的求解方法二分法和Newton法及改进的Newton法。 前言：（目的和意义） 掌握二分法与Newton法的基本原理和应用。 数学原理： 对于一个非线性方程的数值解法很多。在此介绍两种最常见的方法：二分法和Newton法。 对于二分法，其数学实质就是说对于给定的待求解的方程f(x)，其在[a,b]上连续，f(a)f(b)<0，且f(x)在[a,b]内仅有一个实根x*，取区间中点c，若，则c恰为其根，否则根据f(a)f(c)<0是否成立判断根在区间[a,c]和[c,b]中的哪一个，从而得出新区间，仍称为[a,b]。重复运行计算，直至满足精度为止。这就是二分法的计算思想。
Newton法通常预先要给出一个猜测初值x0，然后根据其迭代公式 产生逼近解x*的迭代数列{x k}，这就是Newton法的思想。当x0接近x*时收敛很快，但是当x0选择不好时，可能会发散，因此初值的选取很重要。另外，若将该迭代公式改进为 其中r为要求的方程的根的重数，这就是改进的Newton法，当求解已知重数的方程的根时，在同种条件下其收敛速度要比Newton法快的多。 程序设计： 本实验采用Matlab的M文件编写。其中待求解的方程写成function的方式，如下 function y=f(x); y=-x*x-sin(x); 写成如上形式即可，下面给出主程序。 二分法源程序： clear %%%给定求解区间 b=1.5; a=0;
%%%误差 R=1; k=0;%迭代次数初值 while (R>5e-6) ; c=(a+b)/2; if f12(a)*f12(c)>0; a=c; else b=c; end R=b-a;%求出误差 k=k+1; end x=c%给出解 Newton法及改进的Newton法源程序：clear %%%% 输入函数 f=input('请输入需要求解函数>>','s') %%%求解f(x)的导数 df=diff(f);
学生实验报告实验课程名称 开课实验室 学院年级专业班 学生姓名学号 开课时间至学年学期
if(A(m,k)~=0) if(m~=k) A([k m],:)=A([m k],:); %换行 end A(k+1:n, k:c)=A(k+1:n, k:c)-(A(k+1:n,k)/ A(k,k))*A(k, k:c); %消去end end x=zeros(length(b),1); %回代求解 x(n)=A(n,c)/A(n,n); for k=n-1:-1:1 x(k)=(A(k,c)-A(k,k+1:n)*x(k+1:n))/A(k,k); end y=x; format short;%设置为默认格式显示,显示5位 （2）建立MATLAB界面 利用MA TLAB的GUI建立如下界面求解线性方程组： 详见程序。 五、计算实例、数据、结果、分析 下面我们对以上的结果进行测试，求解：
? ? ? ? ? ? ? ? ? ? ? ? - = ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? - - - - - - 7 2 5 10 13 9 14 4 4 3 2 1 13 12 4 3 3 10 2 4 3 2 1 x x x x 输入数据后点击和，得到如下结果： 更改以上数据进行测试，求解如下方程组： 1 2 3 4 43211 34321 23431 12341 x x x x ?? ???? ?? ???? ?? ???? = ?? ???? - ?? ???? - ???? ?? 得到如下结果：
数值分析实验报告模板 篇一：数值分析实验报告(一)(完整) 数值分析实验报告 1 2 3 4 5 篇二：数值分析实验报告 实验报告一 题目：非线性方程求解 摘要：非线性方程的解析解通常很难给出，因此线性方程的数值解法就尤为重要。本实验采用两种常见的求解方法二分法和Newton法及改进的Newton法。利用二分法求解给定非线性方程的根，在给定的范围内，假设f(x,y)在[a,b]上连续，f(a)xf(b) 直接影响迭代的次数甚至迭代的收敛与发散。即若x0 偏离所求根较远，Newton法可能发散的结论。并且本实验中还利用利用改进的Newton法求解同样的方程，且将结果与Newton法的结果比较分析。 前言：（目的和意义） 掌握二分法与Newton法的基本原理和应用。掌握二分法的原理，验证二分法，在选对有根区间的前提下，必是收
敛，但精度不够。熟悉Matlab语言编程，学习编程要点。体会Newton使用时的优点，和局部收敛性，而在初值选取不当时，会发散。 数学原理： 对于一个非线性方程的数值解法很多。在此介绍两种最常见的方法：二分法和Newton法。 对于二分法，其数学实质就是说对于给定的待求解的方程f(x)，其在[a,b]上连续，f(a)f(b) Newton法通常预先要给出一个猜测初值x0，然后根据其迭代公式xk?1?xk?f(xk) f'(xk) 产生逼近解x*的迭代数列{xk}，这就是Newton法的思想。当x0接近x*时收敛很快，但是当x0选择不好时，可能会发散，因此初值的选取很重要。另外，若将该迭代公式改进为 xk?1?xk?rf(xk) 'f(xk) 其中r为要求的方程的根的重数，这就是改进的Newton 法，当求解已知重数的方程的根时，在同种条件下其收敛速度要比Newton法快的多。 程序设计： 本实验采用Matlab的M文件编写。其中待求解的方程写成function的方式，如下 function y=f(x);
---------------------------------------------------------------最新资料推荐------------------------------------------------------ 数值分析上机实验最小二乘法 数值分析实验报告五最小二乘法一、数值分析实验报告五最 小二乘法一、题目设有如下数据题目设有如下数据 xj -3 -2 -1 0 1 2 3 -3 -2 -1 0 1 2 3 ( )jf x -1.76 0.42 1.2 1.34 1.43 2.25 4.38 -1.76 0.42 1.2 1.34 1.43 2.25 4.38 用三次多项式拟 合这组数据，并绘出图形。 二、用三次多项式拟合这组数据，并绘出图形。 二、方法最小二t(f,[xx(1),xx(n)]) 四、结果 save and run 之后： 请输入插值节点 as [x1,x2...] [-3 -2 -1 0 1 2 3] 请输入插 值节点处对应的函数值 as [f1,f2...] [-1.76 0.42 1.2 1.34 1.43 2.25 4.38] 请输入要求的插值次数m =3 f = 133/100+121469856021/35184372088832*x-8042142191733/450359 结果 save and run 之后： 请输入插值节点 as [x1,x2...] [-3 -2 -1 0 1 2 3] 请输入插 值节点处对应的函数值 as [f1,f2...] [-1.76 0.42 1.2 1.34 1.43 2.25 4.38] 请输入要求的插值次数m =3 f = 133/100+121469856021/35184372088832*x-8042142191733/4503599 627370496*x+1020815915537309/9007199254740992*x 9627370496*x+1020815915537309/9007199254740992*x五、拓展： 1 / 2
实验报告 实验项目名称 估计水塔的水流量 实验项目类型 验证 演示 综合 设计 其他 指导教师 成绩 一、实验目的 （1）学会对实际问题的分析方法 （2）学会利用所学的知识解决实际问题 （3）设计出相应的算法，编制相应的应用程序 二、实验内容 某居民区，其自来水是有一个圆柱形水塔提供，水塔高12.2m ，塔的直径为17.4m ，水塔是由水泵根据水塔中的水位自动加水，一般水泵每天工作两次。按照设计，当水塔中的水位降低至最低水位，约8.2m 时，水泵自动启动加水。当水位升至最高水位，约10.8m 时，水泵停止工作。 下表给出了某一天的测量记录，测量了28个时刻的数据，但由于水泵正向水塔供水，由3个时刻无法测量到水位（表中为—）。 试建立数学模型，计算居民的用水速度和日总用水量。 三、实验原理、方法（算法）、步骤 时刻 0 0.921 1.843 2.949 3.871 4.978 5.900 水位 9.677 9.479 9.308 9.125 8.982 8.814. 8.686 时刻 7.006 7.928 8.967 9.981 10.925 10.945 12.032 水位 8.525 8.388 8.220 — — 10.820 10.500 时刻 12.954 13.875 14.982 15.903 1 6.826 1 7.931 19.037 水位 10.210 9.936 9.653 9.409 9.180 8.921 8.662 时刻 1 9.959 20.839 22.015 22.958 23.880 24.986 25.908 水位 8.433 8.220 — 10.820 10.591 10.354 10.180
序言 数值分析是计算数学的范畴，有时也称它为计算数学、计算方法、数值方法等，其研究对象是各种数学问题的数值方法的设计、分析及其有关的数学理论和具体实现的一门学科，它是一个数学分支。是科学与工程计算（科学计算）的理论支持。许多科学与工程实际问题（核武器的研制、导弹的发射、气象预报）的解决都离不开科学计算。目前，试验、理论、计算已成为人类进行科学活动的三大方法。 数值分析是计算数学的一个主要部分,计算数学是数学科学的一个分支,它研究用计算机求解各种数学问题的数值计算方法及其理论与软件实现。现在面向数值分析问题的计算机软件有：C,C++,MATLAB,Python,Fortran等。 MATLAB是matrix laboratory的英文缩写，它是由美国Mathwork公司于1967年推出的适合用于不同规格计算机和各种操纵系统的数学软件包，现已发展成为一种功能强大的计算机语言，特别适合用于科学和工程计算。目前，MATLAB应用非常广泛，主要用于算法开发、数据可视化、数值计算和数据分析等，除具备卓越的数值计算能力外，它还提供了专业水平的符号计算，文字处理，可视化建模仿真和实时控制等功能。 本实验报告使用了MATLAB软件。对不动点迭代，函数逼近（lagrange插值，三次样条插值，最小二乘拟合），追赶法求解矩阵的解，4RungeKutta方法求解，欧拉法及改进欧拉法等算法做了简单的计算模拟实践。并比较了各种算法的优劣性，得到了对数值分析这们学科良好的理解，对以后的科研数值分析能力有了极大的提高。
目录 序言 (1) 问题一非线性方程数值解法 (3) 1.1 计算题目 (3) 1.2 迭代法分析 (3) 1.3计算结果分析及结论 (4) 问题二追赶法解三对角矩阵 (5) 2.1 问题 (5) 2.2 问题分析（追赶法） (6) 2.3 计算结果 (7) 问题三函数拟合 (7) 3.1 计算题目 (7) 3.2 题目分析 (7) 3.3 结果比较 (12) 问题四欧拉法解微分方程 (14) 4.1 计算题目 (14) 4.2.1 方程的准确解 (14) 4.2.2 Euler方法求解 (14) 4.2.3改进欧拉方法 (16) 问题五四阶龙格-库塔计算常微分方程初值问题 (17) 5.1 计算题目 (17) 5.2 四阶龙格-库塔方法分析 (18) 5.3 程序流程图 (18) 5.4 标准四阶Runge-Kutta法Matlab实现 (19) 5.5 计算结果及比较 (20) 问题六舍入误差观察 (22) 6.1 计算题目 (22) 6.2 计算结果 (22) 6.3 结论 (23) 7 总结 (24) 附录
实验报告 课程名称数值分析 实验项目名称解线性方程组 实验类型上机实验学时 4 班级20111131 学号2011113130 姓名张振指导教师沈艳 实验室名称理学楼407 实验时间2013.12.9 实验成绩预习部分 实验过程 表现 实验报告 部分 总成绩 教师签字日期 哈尔滨工程大学教务处制
实验四 解线性方程组 一．解线性方程组的基本思想 1．直接三角分解法： 将系数矩阵A 转变成等价两个矩阵L 和U 的乘积 ，其中L 和U 分别是下三角和上三角矩阵。当A 的所有顺序主子式都不为0时，矩阵A 可以分解为A=LU ，且分解唯一。其中L 是单位下三角矩阵，U 是上三角矩阵。 2．平方根法： 如果矩阵A 为n 阶对称正定矩阵，则存在一个对角元素为正数的下三角实矩阵L ，使得：A=LL^T 。当限定L 的对角元素为正时，这种分解是唯一的，称为平方根法（Cholesky ）分解。 3．追赶法： 设系数矩阵为三对角矩阵 112223311 10000000000000000 n n n n n b c a b c a b A a b c a b ---?? ? ? ?= ? ? ? ? ?? ? 则方程组Ax=f 称为三对角方程组。 设矩阵A 非奇异，A 有Crout 分解A=LU ，其中L 为下三角矩阵，U 为单位上三角矩阵，记 1122 23311000010 000001000 0000100,0000000000 0001n n n n b L U γαβγββγβ--???? ? ? ? ? ? ? ?== ? ? ? ? ? ? ? ? ? ??? ? ? ? 可先依次求出L ，U 中的元素后，令Ux=y ，先求解下三角方程组Ly=f 得出y ，再求解上三 角方程组Ux=y 。 4．雅克比迭代法： 首先将方程组中的系数矩阵A 分解成三部分，即：A = L+D+U ，如图1所示，其中D 为对角阵，L 为下三角矩阵，U 为上三角矩阵。
数值分析实验(2)
实验二 插值法 P50 专业班级：信计131班 姓名：段雨博 学号：2013014907 一、实验目的 1、熟悉MATLAB 编程； 2、学习插值方法及程序设计算法。 二、实验题目 1、已知函数在下列各点的值为 试用4次牛顿插值多项式()4P x 及三次样条函数()S x （自然边界条件）对数据进行插值用图给出(){},,0.20.08,0,1,11,10i i i x y x i i =+=，()4P x 及()S x 。 2、在区间[]1,1-上分别取10,20n =用两组等距节点对龙格函数()2 1125f x x = +作多项式插值及三次样条插值，对每个n 值，分别画出插值函数及()f x 的图形。 3、下列数据点的插值 可以得到平方根函数的近似，在区间[]0,64上作图 （1）用这9个点作8次多项式插值()8L x （2）用三次样条（第一边界条件）程序求()S x 从得到结果看在[]0,64上，哪个插值更精确；在区间[]0,1上，两种插值哪个更精确？ 三、实验原理与理论基础
1、拉格朗日差值公式 )()(111k k k k k k x x x x y y y x L ---+ =++ 点斜式 k k k k k k k k x x x x y x x x x y x L --+--=++++11111)( 两点式 2、n 次插值基函数 ....,2,1,0,)()(0n j y x l y x L i j n k k k j n ===∑= n k x x x x x x x x x x x x x l n k n k k k k k ,...,1,0,) () (... ) () (... ) () ()(1100=------= -- 3、牛顿插值多项式 ...))(](,,[)](,[)()(102100100+--+++=x x x x x x x f x x x x f x f x P n ))...(](,...,[100---+n n x x x x x x f )(],...,,[)()()(10x x x x f x P x f x R n n n n +=-=ω 4、三次样条函数 若函数],,[)(2b a C x S ∈且在每个小区间],[1+j j x x 上是三次多项式，其中， b x x x a n =<<<=...10是给定节点，则称)(x S 是节点n x x x ,...,,10上的三次样条函数。若在节点j x 上给定函数值),,...,2,1,0)((n j x f y j i ==并成立,,...,2,1,0,)(n j y x S i j ==则称)(x S 为三次样条插值函数。 5、三次样条函数的边界条件 （1）0)()(''''''00''====n n f x S f x S （2）'''00')(,)(n n f x S f x S == 四、实验内容 1、M 文件： function [p]=Newton_Polyfit(X,Y) format long g r=size(X); n=r(2); M=ones(n,n); M(:,1)=Y'; for i=2:n