数值分析上机第二次作业

function [polyout,delta,appval,C]=Legendre_public(fun,a,b,series,xx)

% 该函数为使用Legendre正交多项式求取最佳平方逼近函数(Least Square Approximation) % fun:逼近函数句柄; a,b:逼近函数区间; series:Legendre逼近函数的级数,scalar or vector

% xx:代求函数点; appval:代求函数点的值; polyout:最终求出多项式形式(符号解); delta:逼近函数误差;

% C:对应的Legendre正交多项式系数(降幂排列)

%% 判断输入变量,重建series参数

if numel(series)==1 % 判断输入级数是否为标量

polycoe=zeros(1,series+1); % 仅仅一个值的情况下,polycoe赋初值

C=zeros(1,series+1); % Legendre正交多项式系数赋初值

series1=series:-1:0; % 将单个值转化为级数向量，函数中调用

else polycoe=zeros(1,max(series)+1); % 向量条件下,polycoe赋初值

C=zeros(1,length(series)); % 向量条件下，Legendre正交多项式系数赋初值

series1=sort(series,'descend'); % 保证输入的向量级数成降序排列

end

%% 计算正交多项式系数Cj & 误差余项delta

delta_handle=@(t)delta_integral(t,fun,a,b); % trap_integral函数中第一个参数

delta=trap_integral(delta_handle,-1,1,1e-3); % 误差余项delta赋初值

for jj=1:size(C,2)

C_handle=@(t)f_integral(t,fun,a,b,series1(jj));

C(jj)=(2*series1(jj)+1)/2*(trap_integral(C_handle,-1,1)); % 计算级数,Legendre正交多项式级数降幂排列

delta=delta-C(jj)^2*2/(2*series1(jj)+1); % 求取误差余项delta

end

%% 计算polycoe多项式系数(降幂排列)

Coe_cell=cell(length(series1),1); % 创建胞元数组,存储Legendre级数

for ii=1:length(series1) % Legendre级数进行迭代

Coe_cell{ii}=tailassign(Legen_poly(series1(ii)),series1(1)+1); % tailassign函数作用不同级数Legendre多项式

polycoe=C(ii)*Coe_cell{ii}+polycoe; % polycoe向量的值,Legendre多项式Pn(t)表示end

string=[num2str(2/(b-a)),'*','x','+',num2str((b+a)/(a-b))];

polyexp=simplify(subs(poly2sym(polycoe,'t'),'t',string)); % 符号函数作符号变化

%% 计算待求点的拟合函数值

if nargin==5

appval=double(subs(polyexp,xx)); % 利用符号函数的替换计算待求点的数值

else xx=0; appval=double(subs(polyexp,xx)); % 若输入变量不足5个,xx默认赋值

end

polyout=vpa(polyexp,5);

%% 内嵌函数1:构造Legendre正交多项式

function Legen_coe=Legen_poly(ser)

% 该函数用于构造Lengendre正交多项式, ser:为正交多项式级数,标量; Lengen_coe:返回多项式系数(降幂排列),向量

if ser==0

Legen_coe=1; % 当Legendre多项式级数为0时的Legendre多项式else

if ser==1

Legen_coe=[1,0]; % 当Legendre多项式级数为1时的Legendre多项式系数else % 当Legendre多项式级数为1时的Legendre多项式系数

Legen_coe=(2*ser-1)/ser*[Legen_poly(ser-1),0]-(ser-1)/ser*[0,0,Legen_poly(ser-2)];

end

end

end

%% 内嵌函数2:构造求取Cj系数需要求积函数

function fpfun=f_integral(t,fhandle,a,b,n)

% 该函数用于构建f(x)*Pj(x)在[-1,1]区间上的求积公式，用于积分求取时的句柄;函数还涉及积分区间变换

% t:区间变化后的值; fhandle:被逼近原函数句柄; a,b:原函数积分区间; n:Legendre多项式Pj(x)级数,scalar;

% fpfun:最终函数表达式

Legpoly=Legen_poly(n);

x0=(b-a)/2*t+(b+a)/2;

fpfun=polyval(Legpoly,t).*feval(fhandle,x0); % 构建f(x)*Pj(t)的积分函数end

%% 内嵌函数3:构造求取误差余项所需要的函数

function ffun=delta_integral(t,fhandle,a,b)

% 该函数用于构建f(x)*f(x)在[-1,1]区间上的求积公式，用于积分求取时的句柄;函数还涉及积分区间变换

% t:区间变化后的值; fhandle:被逼近原函数句柄; a,b:原函数积分区间;

% ffun:最终函数表达式

x0=(b-a)/2*t+(b+a)/2; % x区间从[a,b]变化到t[-1,1]

ffun=feval(fhandle,x0).*feval(fhandle,x0);

end

%% 内嵌函数4：构建赋值函数,向量直接求取多项式系数(降幂排列)

function Vec=tailassign(Mo,len)

% 该函数用于构建矩阵延后赋值函数,Vec向量：Mo向量中的元素首项向前填充0直到满足len长度为止

% Vec:构造降幂排列的多项式系数; Mo:初始向量,Mo向量长度

Vec=zeros(1,len);

account=1;

if len

error('尾部赋值函数,输出向量长度要不小于输入向量长度')

end

for ki=len-length(Mo)+1:len

Vec(ki)=Mo(account); % 将Mo向量向后充填在Vec向量中

account=account+1;

end

end

end

function Tn2=trap_integral(handle,a,b,tol)

% 该函数：变步长复化梯形公式求取函数积分

% handle:函数句柄; a,b:积分函数上下限; tol：函数允许绝对误差

% Tn2:满足误差条件下的数值解

%% 输入条件判断& 设置函数使用参数

if nargin<4 || isempty(tol)

tol=1e-3;

end

iteration_max=50; % 设置最大迭代次数

%% 循环迭代

iteration=1; % 设置初始迭代次数

T0=trap_varstep(handle,a,b,1,iteration);

T1=trap_varstep(handle,a,b,T0,iteration+1);

iteration=iteration+1;

while abs(T1-T0)>tol && iteration

T0=T1;

T1=trap_varstep(handle,a,b,T0,iteration+1);

sprintf('第%d 次迭代，迭代误差: %d',iteration,abs(T1-T0))

iteration=iteration+1;

end

%% 结果判断

if iteration

sprintf('迭代收敛, 最终收敛误差为：%d',abs(T1-T0))

Tn2=T1;

else warning('迭代不收敛, 最终收敛误差为：%d',abs(T1-T0))

Tn2=T1;

end

%% 变步长复化梯形公式求取积分

function T2=trap_varstep(handle,a,b,T0,iter)

% 该函数用于计算变步长梯形公式

% handle:函数句柄; a,b:积分上下限; T0:变步长递推公式下的原积分值;iter:积分上下限iter等分

% T2:满足误差限度下的积分值

if iter==1 || iter==0

f0=feval(handle,a); f1=feval(handle,b);

T2=(f0+f1)*(b-a)/2;

else

n=2^(iter-1);

h=(b-a)/n;

xin=a:h:b;

Hno=0;

for jj=2:length(xin)-1

ff=feval(handle,xin(jj));

Hno=Hno+ff;

T2=1/2*(T0+Hno*h);

if ~mod(jj,1000000)

sprintf('计算复化梯形中节点和，当前迭代%d 次, 当前中间节点和为: %d',jj,Hno)

end

end

end

end

end

a=0;b=4;

f=@(x)(cos(x).*exp(x));

[coe,delta,app,C]=Legendre_public(f,a,b,3);

xx=linspace(0,4,1e3+1);

yy=subs(coe,'x',xx);

yval=double(yy);

[coe1,delta1,app1,C1]=Legendre_public(f,a,b,4);

yy1=subs(coe1,'x',xx);

yval2=double(yy1);

yreal=feval(f,xx);

%————————————做图————————————————————figure

plot(xx,yval2,'bh','Markersize',5)

hold on

plot(xx,yval2,'rs','Markersize',2)

plot(xx,yreal,'g--','Linewidth',3)

grid on

xlabel('X轴')

ylabel('Y轴')

legend('Legendre正交多项式:级数4','Legendre正交多项式:级数3','原函数')

hold off

数值分析上机作业

昆明理工大学工科研究生《数值分析》上机实验 学院：材料科学与工程学院 专业：材料物理与化学 学号：2011230024 姓名: 郑录 任课教师：胡杰

P277-E1 1.已知矩阵A= 10787 7565 86109 75910 ?? ?? ?? ?? ?? ??,B= 23456 44567 03678 00289 00010 ?? ?? ?? ?? ?? ?? ?? ?? ,错误！未找到引用源。 = 11/21/31/41/51/6 1/21/31/41/51/61/7 1/31/41/51/61/71/8 1/41/51/61/71/81/9 1/51/61/71/81/91/10 1/61/71/81/91/101/11?????????????????? （1）用MA TLAB函数“eig”求矩阵全部特征值。 （2）用基本QR算法求全部特征值（可用MA TLAB函数“qr”实现矩阵的QR分解）。解：MA TLAB程序如下： 求矩阵A的特征值： clear; A=[10 7 8 7;7 5 6 5;8 6 10 9;7 5 9 10]; E=eig(A) 输出结果： 求矩阵B的特征值： clear; B=[2 3 4 5 6;4 4 5 6 7;0 3 6 7 8;0 0 2 8 9;0 0 0 1 0]; E=eig(B) 输出结果：

求矩阵错误！未找到引用源。的特征值： clear; 错误！未找到引用源。=[1 1/2 1/3 1/4 1/5 1/6; 1/2 1/3 1/4 1/5 1/6 1/7; 1/3 1/4 1/5 1/6 1/7 1/8; 1/4 1/5 1/6 1/7 1/8 1/9;1/5 1/6 1/7 1/8 1/9 1/10; 1/6 1/7 1/8 1/9 1/10 1/11]; E=eig(错误！未找到引用源。) 输出结果： （2）A= 10 7877565861097 5 9 10 第一步：A0=hess(A);[Q0,R0]=qr(A0);A1=R0*Q0 返回得到： 第二部：[Q1,R1]=qr(A1);A2=R1*Q1

数值分析上机作业

数值分析上机实验报告 选题：曲线拟合的最小二乘法 指导老师： 专业： 学号： 姓名：

课题八曲线拟合的最小二乘法 一、问题提出 从随机的数据中找出其规律性，给出其近似表达式的问题，在生产实践和科学实验中大量存在，通常利用数据的最小二乘法求得拟合曲线。 在某冶炼过程中，根据统计数据的含碳量与时间关系，试求含碳量y 与时间t 的拟合曲线。 二、要求 1、用最小二乘法进行曲线拟合； 2、近似解析表达式为()33221t a t a t a t ++=?； 3、打印出拟合函数()t ?，并打印出()j t ?与()j t y 的误差，12,,2,1 =j ； 4、另外选取一个近似表达式，尝试拟合效果的比较； 5、*绘制出曲线拟合图*。 三、目的和意义 1、掌握曲线拟合的最小二乘法； 2、最小二乘法亦可用于解超定线代数方程组； 3、探索拟合函数的选择与拟合精度间的关系。 四、计算公式 对于给定的测量数据(x i ,f i )(i=1,2,…，n )，设函数分布为 ∑==m j j j x a x y 0)()(? 特别的，取)(x j ?为多项式 j j x x =)(? (j=0, 1,…，m )

则根据最小二乘法原理，可以构造泛函 ∑∑==-=n i m j i j j i m x a f a a a H 1 10))((),,,(? 令 0=??k a H (k=0, 1,…，m ) 则可以得到法方程 ???? ??????? ?=????????????????????????),(),(),(),(),(),(),(),(),(),(),(),(1010101111000100m m m m m m m m f f f a a a ????????????????????? 求该解方程组，则可以得到解m a a a ,,,10 ，因此可得到数据的最小二乘解 ∑=≈m j j j x a x f 0)()(? 曲线拟合：实际工作中，变量间未必都有线性关系，如服药后血药浓度与时间的关系；疾病疗效与疗程长短的关系；毒物剂量与致死率的关系等常呈曲线关系。曲线拟合是指选择适当的曲线类型来拟合观测数据，并用拟合的曲线方程分析两变量间的关系。 五、结构程序设计 在程序结构方面主要是按照顺序结构进行设计，在进行曲线的拟合时，为了进行比较，在程序设计中，直接调用了最小二乘法的拟合函数polyfit ，并且依次调用了plot 、figure 、hold on 函数进行图象的绘制，最后调用了一个绝对值函数abs 用于计算拟合函数与原有数据的误差，进行拟合效果的比较。

东南大学数值分析上机题答案

数值分析上机题 第一章 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 运行结果：

东南大学数值分析上机作业汇总

东南大学数值分析上机作业 汇总 -标准化文件发布号：（9456-EUATWK-MWUB-WUNN-INNUL-DDQTY-KII

数值分析上机报告 院系： 学号： 姓名：

目录 作业1、舍入误差与有效数 (1) 1、函数文件cxdd.m (1) 2、函数文件cddx.m (1) 3、两种方法有效位数对比 (1) 4、心得 (2) 作业2、Newton迭代法 (2) 1、通用程序函数文件 (3) 2、局部收敛性 (4) （1）最大δ值文件 (4) （2）验证局部收敛性 (4) 3、心得 (6) 作业3、列主元素Gauss消去法 (7) 1、列主元Gauss消去法的通用程序 (7) 2、解题中线性方程组 (7) 3、心得 (9) 作业4、三次样条插值函数 (10) 1、第一型三次样条插值函数通用程序： (10) 2、数据输入及计算结果 (12)

作业1、舍入误差与有效数 设∑ =-=N j N j S 2 2 11 ，其精确值为?? ? ??---1112321N N . （1）编制按从小到大的顺序1 1 131121222-? ??+-+-=N S N ，计算N S 的通用程序； （2）编制按从大到小的顺序()1 21 11111222-???+--+-=N N S N ，计算N S 的通用程序； （3）按两种顺序分别计算642101010,,S S S ,并指出有效位数； （4）通过本上机你明白了什么？ 程序： 1、函数文件cxdd.m function S=cxdd(N) S=0; i=2.0; while (i<=N) S=S+1.0/(i*i-1); i=i+1; end script 运行结果（省略>>）： S=cxdd(80) S= 0.737577 2、函数文件cddx.m function S=cddx (N) S=0; for i=N:-1:2 S=S+1/(i*i-1); end script 运行结果（省略>>）： S=cddx(80) S= 0.737577 3、两种方法有效位数对比

数值分析上机题目详解

第一章 一、题目 设∑ =-= N N j S 2 j 2 1 1，其精确值为)11 123(21+--N N 。 1) 编制按从大到小的顺序1 1 13112122 2-+??+-+-=N S N ，计算S N 的通用程序。 2) 编制按从小到大的顺序1 21 1)1(111222-+ ??+--+-= N N S N ，计算S N 的通用程序。 3) 按两种顺序分别计算64210,10,10S S S ,并指出有效位数。（编制程序时用单精度） 4) 通过本次上机题，你明白了什么？ 二、通用程序 N=input('Please Input an N (N>1):'); AccurateValue=single((0-1/(N+1)-1/N+3/2)/2); Sn1=single(0); for a=2:N; Sn1=Sn1+1/(a^2-1); end Sn2=single(0); for a=2:N; Sn2=Sn2+1/((N-a+2)^2-1); end fprintf('The value of Sn (N=%d)\n',N); fprintf('Accurate Calculation %f\n',AccurateValue); fprintf('Caculate from large to small %f\n',Sn1); fprintf('Caculate from small to large %f\n',Sn2); disp('____________________________________________________')

三、结果 从结果可以看出有效位数是6位。 感想：可以得出，算法对误差的传播有一定的影响，在计算时选一种好的算法可以使结果更为精确。从以上的结果可以看到从大到小的顺序导致大数吃小数的现象，容易产生较大的误差，求和运算从小数到大数所得到的结果才比较准确。

Matlab作业3(数值分析)答案

Matlab作业3（数值分析） 机电工程学院（院、系）专业班组 学号姓名实验日期教师评定 1.计算多项式乘法(x2+2x+2)(x2+5x+4)。 答： 2. （1）将(x-6)(x-3)(x-8)展开为系数多项式的形式。（2）求解在x=8时多项 式(x-1)(x-2) (x-3)(x-4)的值。 答：（1） （2）

3. y=sin(x)，x从0到2π，?x=0.02π，求y的最大值、最小值、均值和标准差。 4.设ｘ＝[0.00.30.8 1.1 1.6 2.3]'，y=[0.500.82 1.14 1.25 1.35 1.40]'，试求二次多项式拟合系数，并据此计算x1=[0.9 1.2]时对应的y1。解：x=[0.0 0.3 0.8 1.1 1.6 2.3]'; %输入变量数据x y=[0.50 0.82 1.14 1.25 1.35 1.40]'; %输入变量数据y p=polyfit(x,y,2) %对x,y用二次多项式拟合,得到系数p x1=[0.9 1.2]; %输入点x1 y1=polyval(p,x1) %估计x1处对应的y1 p = -0.2387 0.9191 0.5318 y1 = a) 1.2909

5.实验数据处理：已知某压力传感器的测试数据如下表 p为压力值，u为电压值，试用多项式 d cp bp ap p u+ + + =2 3 ) ( 来拟 合其特性函数，求出a,b,c,d，并把拟合曲线和各个测试数据点画在同一幅图上。解： >> p=[0.0,1.1,2.1,2.8,4.2,5.0,6.1,6.9,8.1,9.0,9.9]; u=[10,11,13,14,17,18,22,24,29,34,39]; x=polyfit(p,u,3) %得多项式系数 t=linspace(0,10,100); y=polyval(x,t); %求多项式得值 plot(p,u,'*',t,y,'r') %画拟和曲线 x = 0.0195 -0.0412 1.4469 9.8267

2014级硕士研究生数值分析上机实习报告

2014级硕士研究生数值分析上机实习（第一次） 姓名：学号：学院： 实习题目：分别用二分法和Newton迭代法求方程x3■ 2x210x-20=0的根.实习目的：掌握两种解法，体会两种解法的收敛速度. 实习要求：用C程序语言编程上机进行计算，精确到8位有效数字. 报告内容： 1.确定实根的个数以及所在区间 2.将最后两次计算结果填入下表（保留8位数字）： 3.实习过程中遇到哪些问题？如何解决？有何心得体会?

4.两种解法的计算程序（此页写不下时可以加页）：

2014级硕士研究生数值分析上机实习（第二次）姓名：学号：学院： 实习题目：计算8阶三对角矩阵A=tridiag（0.235, 1.274, 0.235）的行列式.实习目的：掌握计算行列式的方法. 实习要求：首先选择一种算法，然后用C程序语言编程上机进行计算.报告内容： 1.简单描述所采用的算法: 2?计算结果： A 3.实习过程中遇到哪些问题？如何解决？有何心得体会?

4.写出C语言计算程序（此页写不下时可以加页）：

2014级硕士研究生数值分析上机实习（第三次） 姓名：学号：学院： 分别用Jacobi迭代法和Gauss-Seidel迭代法求解线性方程组实习题目： 2lx + 9.8y+ 3.4z= 6.7 <2.7x + 1.8y+ 7.2z= 2.4 8.6x + 1.5y + 3.4z = 1.9 实习目的：感受两种迭代法的收敛速度. 首先构造收敛的Jacobi迭代法和Gauss-Seidel迭代法，然后用实习要求： C程序语言编程上机进行求解，初始值均取为0,精确到4位小 数. 报告内容： 1.写出收敛的Jacobi迭代法和Gauss-Seidel迭代法:

北航数值分析大作业第二题精解

目标：使用带双步位移的QR 分解法求矩阵10*10[]ij A a =的全部特征值，并对其中的每一个实特征值求相应的特征向量。已知：sin(0.50.2)() 1.5cos( 1.2)(){i j i j ij i j i j a +≠+== (i,j=1,2, (10) 算法： 以上是程序运作的逻辑，其中具体的函数的算法，大部分都是数值分析课本上的逻辑，在这里特别写出矩阵A 的实特征值对应的一个特征向量的求法： ()[]()() []()[]()111111I 00000 i n n n B A I gause i n Q A I u Bu u λλ-?-?-=-?-?? ?-=????→=??????→= ?? ? 选主元的消元 检查知无重特征值 由于=0i A I λ- ，因此在经过选主元的高斯消元以后，i A I λ- 即B 的最后一行必然为零，左上方变 为n-1阶单位矩阵[]()()11I n n -?-，右上方变为n-1阶向量[]()11n Q ?-，然后令n u 1=-，则 ()1,2,,1j j u Q j n ==???-。

这样即求出所有A所有实特征值对应的一个特征向量。 #include #include #include #define N 10 #define E 1.0e-12 #define MAX 10000 //以下是符号函数 double sgn(double a) { double z; if(a>E) z=1; else z=-1; return z; } //以下是矩阵的拟三角分解 void nishangsanjiaodiv(double A[N][N]) { int i,j,k; int m=0; double d,c,h,t; double u[N],p[N],q[N],w[N]; for(i=0;i

《数值分析》上机实验报告

数值分析上机实验报告 精选

《数值分析》上机实验报告 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 ==- ==∈≤-≠>+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 #include main() {double x2,f,f1; double x1=1.9; //取初值为1.9 do {x2=x1; f=pow(x2,7)-28*pow(x2,4)+14; f1=7*pow(x2,6)-4*28*pow(x2,3); x1=x2-f/f1;} while(fabs(x1-x2)>=0.00001||x1<0.1); //限制循环次数printf("计算结果：x=%f\n",x1);} 1.3 运行结果： 1.4 MATLAB上机程序 function y=Newton(f,df,x0,eps,M) d=0; for k=1:M if feval(df,x0)==0 d=2;break else x1=x0-feval(f,x0)/feval(df,x0); end e=abs(x1-x0); x0=x1; if e<=eps&&abs(feval(f,x1))<=eps d=1;break end end

东南大学-数值分析上机题作业-MATLAB版

2015.1.9 上机作业题报告 JONMMX 2000

1.Chapter 1 1.1题目 设S N =∑1j 2?1 N j=2 ，其精确值为 )1 1 123(21+--N N 。 （1）编制按从大到小的顺序1 1 131121222-+ ??+-+-=N S N ，计算S N 的通用程序。 （2）编制按从小到大的顺序1 21 1)1(111222-+ ??+--+-= N N S N ，计算S N 的通用程序。 （3）按两种顺序分别计算64210,10,10S S S ,并指出有效位数。（编制程序时用单精度） （4）通过本次上机题，你明白了什么？ 1.2程序 1.3运行结果

1.4结果分析 按从大到小的顺序，有效位数分别为：6，4，3。 按从小到大的顺序，有效位数分别为：5，6，6。 可以看出，不同的算法造成的误差限是不同的，好的算法可以让结果更加精确。当采用从大到小的顺序累加的算法时，误差限随着N 的增大而增大，可见在累加的过程中，误差在放大，造成结果的误差较大。因此，采取从小到大的顺序累加得到的结果更加精确。 2.Chapter 2 2.1题目 （1）给定初值0x 及容许误差ε，编制牛顿法解方程f(x)=0的通用程序。 （2）给定方程03 )(3 =-=x x x f ,易知其有三个根3,0,3321= *=*-=*x x x ○1由牛顿方法的局部收敛性可知存在,0>δ当),(0δδ+-∈x 时，Newton 迭代序列收敛于根x2*。试确定尽可能大的δ。 ○2试取若干初始值，观察当),1(),1,(),,(),,1(),1,(0+∞+-----∞∈δδδδx 时Newton 序列的收敛性以及收敛于哪一个根。 （3）通过本上机题，你明白了什么？ 2.2程序

数值分析第二次大作业

《数值分析》计算实习报告 第二题 院系：机械工程及自动化学院 学号： 姓名： 2017年11月

一、题目要求 试求矩阵A =[a ij ]10×10的全部特征值，并对其中的每一个实特征值求相应的特征向量，已知 a ij ={ sin (0.5i +0.2j ) i ≠j 1.52cos (i +1.2j ) i =j (i,j =1,2, (10) 说明： 1.用带双步位移的QR 方法求矩阵特征值，要求迭代的精度水平为ε=10?12。 2.打印以下内容： （1）全部源程序； （2）矩阵A 经过拟上三角化后所得的矩阵A (n?1)； （3）对矩阵A (n?1)实行QR 方法迭代结束后所得的矩阵； （4）矩阵A 的全部特征值λi =(R i ,I i ) (i =1,2,?,10)，其中R i =Re(λi )，I i = Im(λi ) 。若λi 是实数，则令I i =0； （5）A 的相应于实特征值的特征向量。 3.采用e 型数输出实型数，并且至少显示12位有效数字。 二、算法设计思路和方案 1. 将矩阵A 拟上三角化得到矩阵A (n?1) 为了减少计算量，一般先利用Householder 矩阵对矩阵A 作相似变换，把A 化为拟上三角矩阵A (n?1)，然后用QR 方法计算A (n?1)的全部特征值，而A (n?1)的特征值就是A 的特征值。具体算法如下： 记(1)A A =，()r A 的第r 列至第n 列的元素为(r)(1,2, ,;,1,,)ij a i n j r r n ==+。 对于1,2,,2r n =-执行 （1）若() (2,3,,)r ir a i r r n =++全为零，则令(1)()r r A A +=，转(5)；否则转（2）。

高等数值分析上机作业

高等数值分析上机作业

目录 上机作业1 (1) 上机作业2 (5) 上机作业3 (10) 上机作业4 (13) 上机作业5 (16) 上机作业6 (19) 上机作业7 (20) 上机作业8 (29)

第8章 函数逼近与曲线拟合 上机作业1： 最佳平方逼近 8-11.设()[]1,1,arcsin -∈=x x x f ， （1） 在{}32,,,1x x x span =φ中求()x f 的最佳平方逼近多项式； （2） 在{})(),(),(),(3210x T x T x T x T span =φ中求()x f 的最佳平方逼近多项式。 解：（1） 基于幂函数的最佳平方逼近 简单原理： 对于],[)(b a C x f ∈及一个线性无关函数组的集合 {}],,[)(,),(),(10b a C x x x span n ?=???φ 若存在,φ∈*S 使得 ()dx x S x f x x S x f x S x f b a x S x S ?-=-=-∈∈* 2)(2 2)(2 2 )]()([min )()(min )()(ρφ φ ，则称()x S *是 ()x f 在子集[]b a ,?φ中的最佳平方逼近函数。 取(),,,1,0,n j x x j j ==?就有{}n x x x span ,,,,12 =φ。对于任意的()φ∈x S ，有()∑==n j j j x a x S 0，()x S 为次数n ≤的多项式。 令)(x f 在},,,1{32x x x span =φ中的最佳平方逼近函数为 φ∈+++=3 *2**1*0*3 2)(x a x a x a a x S 通过求解法方程 ???? ? ? ? ??=??????? ????????? ??),(),(),(),(),(),(),(),(),(),(),(),(),(),(),(),(),(),(),() ,(321032103323130 3322212023121110130201000????????????????????????????????????f f f f a a a a 其中.arcsin )(,)(,)(,)(,1)(332210x x f x x x x x x x =====????

(完整版)数值计算方法上机实习题答案

1． 设?+=1 05dx x x I n n ， （1） 由递推公式n I I n n 1 51+-=-，从0I 的几个近似值出发，计算20I ； 解：易得：0I =ln6-ln5=0.1823, 程序为： I=0.182; for n=1:20 I=(-5)*I+1/n; end I 输出结果为：20I = -3.0666e+010 （2） 粗糙估计20I ，用n I I n n 51 5111+- =--，计算0I ； 因为 0095.05 6 0079.01020 201 020 ≈<<≈??dx x I dx x 所以取0087.0)0095.00079.0(2 1 20=+= I 程序为：I=0.0087; for n=1:20 I=(-1/5)*I+1/(5*n); end I 0I = 0.0083 （3） 分析结果的可靠性及产生此现象的原因(重点分析原因)。 首先分析两种递推式的误差；设第一递推式中开始时的误差为000I I E '-=，递推过程的舍入误差不计。并记n n n I I E '-=，则有01)5(5E E E n n n -==-=-Λ。因为=20E 20020)5(I E >>-，所此递推式不可靠。而在第二种递推式中n n E E E )5 1(5110-==-=Λ，误差在缩小， 所以此递推式是可靠的。出现以上运行结果的主要原因是在构造递推式过程中，考虑误差是否得到控制， 即算法是否数值稳定。 2． 求方程0210=-+x e x 的近似根，要求4 1105-+?<-k k x x ，并比较计算量。 （1） 在[0，1]上用二分法； 程序：a=0;b=1.0; while abs(b-a)>5*1e-4 c=(b+a)/2;

数值分析上机题(matlab版)(东南大学)

数值分析上机题(matlab版)(东南大学)

数值分析上机报告

第一章 一、题目 精确值为)1 1 123(21+--N N 。 1) 编制按从大到小的顺序 1 1 131121222-+??+-+-= N S N ，计算S N 的通用程序。 2) 编制按从小到大的顺序 1 21 1)1(111222-+??+--+-= N N S N ，计算S N 的通用程序。 3) 按两种顺序分别计算6 42 10,10, 10S S S ,并指出有效位 数。（编制程序时用单精度） 4) 通过本次上机题，你明白了什么？ 二、通用程序 clear N=input('Please Input an N (N>1):'); AccurateValue=single((0-1/(N+1)-1/N+3/2)/2); Sn1=single(0); for a=2:N; Sn1=Sn1+1/(a^2-1); end Sn2=single(0); for a=2:N; Sn2=Sn2+1/((N-a+2)^2-1); end fprintf('The value of Sn using different algorithms (N=%d)\n',N); disp('____________________________________________________') fprintf('Accurate Calculation %f\n',AccurateValue); fprintf('Caculate from large to small %f\n',Sn1); fprintf('Caculate from small to large %f\n',Sn2);

数值分析上机作业1-1

数值计算方法上机题目1 1、实验1. 病态问题 实验目的： 算法有“优”与“劣”之分，问题也有“好”和“坏”之别。所谓坏问题就是问题本身的解对数据变化的比较敏感，反之属于好问题。希望读者通过本实验对此有一个初步的体会。 数值分析的大部分研究课题中，如线性代数方程组、矩阵特征值问题、非线性方程及方程组等都存在病态的问题。病态问题要通过研究和构造特殊的算法来解决，当然一般要付出一些代价（如耗用更多的机器时间、占用更多的存储空间等）。 问题提出： 考虑一个高次的代数多项式 ∏=-= ---=20 1)()20)...(2)(1()(k k x x x x x p （E1-1） 显然该多项式的全部根为l ，2，…，20，共计20个，且每个根都是单重的（也称为简 单的）。现考虑该多项式方程的一个扰动 0)(19 =+x x p ε (E1-2) 其中ε是一个非常小的数。这相当于是对（E1-1）中19 x 的系数作一个小的扰动。我们希望比较（E1-1）和（E1-2）根的差别，从而分析方程（E1-1）的解对扰动的敏感性。 实验内容： 为了实现方便，我们先介绍两个 Matlab 函数：“roots ”和“poly ”，输入函数 u ＝roots （a ） 其中若变量a 存储1+n 维的向量，则该函数的输出u 为一个n 维的向量。设a 的元素依次为121,...,,+n a a a ，则输出u 的各分量是多项式方程 0...1121=++++-n n n n a x a x a x a 的全部根，而函数 b=poly(v) 的输出b 是一个n ＋1维变量，它是以n 维变量v 的各分量为根的多项式的系数。可见“roots ”和“Poly ”是两个互逆的运算函数. ve=zeros(1,21); ve(2)=ess; roots(poly(1:20))+ve) 上述简单的Matlab 程序便得到（E1-2）的全部根，程序中的“ess ”即是（E1-2）中的ε。 实验要求： （1）选择充分小的ess ，反复进行上述实验，记录结果的变化并分析它们。如果扰动项的系数ε很小，我们自然感觉（E1-1）和(E1-2)的解应当相差很小。计算中你有什么出乎意料的发现？表明有些解关于如此的扰动敏感性如何？ （2）将方程（E1-2）中的扰动项改成18 x ε或其他形式，实验中又有怎样的现象出现？

北航数值分析第二次大作业--QR分解

《数值分析A》

一、算法设计方案 整个程序主要分为四个函数，主函数，拟上三角化函数，QR分解函数以及使用双步位移求解矩阵特征值、特征向量的函数。因为在最后一个函数中也存在QR分解，所以我没有采用参考书上把矩阵M进行的QR分解与矩阵Ak的迭代合并的方法，而是在该函数中调用了QR分解函数，这样增强了代码的复用性，减少了程序长度；但由于时间关系，对阵中方法的运算速度没有进行深入研究。 1.为了减少QR分解法应用时的迭代次数，首先对给定矩阵进行拟上三角化处理。 2.对经过拟上三角化处理的矩阵进行QR分解。 3.注意到计算特征值与特征向量的过程首先要应用前面两个函数，于是在拟上三角化矩阵的基础上对QR分解函数进行了调用。计算过程中，没有采用goto语句，而是根据流程图采用其他循环方式完成了设计，通过对迭代过程的合并，简化了程序的循环次数，最后在计算特征向量的时候采用了列主元高斯消去法。

二、源程序代码 #include #include #include int i,j,k,l,m; //定义外部变量double d,h,b,c,t,s; double A[10][10],AA[10][10],R[10][10],Q[10][10],RQ[10][10]; double X[10][10],Y[10][10],Qt[10][10],M[10][10]; double U[10],P[10],T[10],W[10],Re[10]={0},Im[10]={0}; double epsilon=1e-12; void main() { void Quasiuppertriangular(double A[][10]); void QRdecomposition(double A[][10]); void DoublestepsQR(double A[][10]); int i,j; for(i=0;i<10;i++) { for(j=0;j<10;j++) { A[i][j]=sin(0.5*(i+1)+0.2*(j+1)); Q[i][j]=0; AA[i][j]=A[i][j]; } A[i][i]=1.5*cos(2.2*(i+1)); AA[i][i]=A[i][i];

数值分析作业

第二章 1. 题目：运用MATLAB编程实现牛顿迭代 2. 实验操作 1、打开MATLAB程序软件。 2、在MATLAB中编辑如下的M程序。 function [p1,err,k,y]=newton(f,df,p0,delta,max) %f 是要求根的方程(f(x)=0); %df 是f(x)的导数； %p0是所给初值,位于x*附近； %delta是给定允许误差； %max是迭代的最大次数； %p1是newton法求得的方程的近似解； %err是p0的误差估计； %k是迭代次数； p0 for k=1:max p1=p0-feval('f',p0)/feval('df',p0); err=abs(p1-p0); p0=p1; k p1 err y=feval('f',p1) if (err> newton('f','df',1.2,10^(-6),20) 3.实验结果

p0 = 1.2000 k =1 p1=1.1030 err=0.0970 y=0.0329 k= 2 p1=1.0524 err=0.0507 y=0.0084 k =3 p1=1.0264 err=0.0260 y=0.0021 k =4 p1=1.0133 err=0.0131 y=5.2963e-004 k =5 p1=1.0066 err=0.0066 y=1.3270e-004 k =6 p1=1.0033 err=0.0033 y=3.3211e-005 k =7 p1=1.0017 err=0.0017 y=8.3074e-006 k =8 p1=1.0008 err=8.3157e-004 y = 2.0774e-006 k =9 p1=1.0004 err=4.1596e-004 y =5.1943e-007 k=10 p1=1.0002 err=2.0802e-004 y= 1.2987e-007 k=11 p1=1.0001 err=1.0402e-004 y =3.2468e-008 k=12 p1=1.0001 err=5.2014e-005 y=8.1170e-009 k=13 p1=1.0000 err=2.6008e-005 y= 2.0293e-009 k=14 p1=1.0000 err=1.3004e-005 y=5.0732e-010 k=15 p1 =1.0000 err=6.5020e-006 y=1.2683e-010 k=16 p1 =1.0000 err=3.2510e-006 y=3.1708e-011 k=17 p1 =1.0000 err=1.6255e-006 y =7.9272e-012 k=18 p1 =1.0000 err =8.1279e-007 y= 1.9820e-012 ans = 1.0000 结果说明：经过18次迭代得到精确解为1，误差为8.1279e-007。

北航数值分析报告大作业第八题

北京航空航天大学 数值分析大作业八 学院名称自动化 专业方向控制工程 学号 学生姓名许阳 教师孙玉泉 日期2014 年11月26 日

一．题目 关于x , y , t , u , v , w 的方程组(A.3) ???? ?? ?=-+++=-+++=-+++=-+++79 .0sin 5.074.3cos 5.007.1cos sin 5.067.2cos 5.0y w v u t x w v u t y w v u t x w v u t (A.3) 以及关于z , t , u 的二维数表（见表A-1）确定了一个二元函数z =f (x , y )。 表A-1 二维数表 t z u 0 0.4 0.8 1.2 1.6 2 0 -0.5 -0.34 0.14 0.94 2.06 3.5 0.2 -0.42 -0.5 -0.26 0.3 1.18 2.38 0.4 -0.18 -0.5 -0.5 -0.18 0.46 1.42 0.6 0.22 -0.34 -0.58 -0.5 -0.1 0.62 0.8 0.78 -0.02 -0.5 -0.66 -0.5 -0.02 1.0 1.5 0.46 -0.26 -0.66 -0.74 -0.5 1. 试用数值方法求出f (x , y ) 在区域}5.15.0,8.00|), {≤≤≤≤=y x y x D （上的近似表达式 ∑∑===k i k j s r rs y x c y x p 00 ),( 要求p (x , y )以最小的k 值达到以下的精度 ∑∑==-≤-=10020 7210)],(),([i j i i i i y x p y x f σ 其中j y i x i i 05.05.0,08.0+==。 2. 计算),(),,(* ***j i j i y x p y x f (i =1,2,…,8 ; j =1,2,…,5) 的值，以观察p (x , y ) 逼 近f (x , y )的效果，其中j y i x j i 2.05.0,1.0**+==。

东南大学数值分析上机题(上)

数值分析上机报告 姓名： 学号： 专业： 2013年10月27日

第一章 舍入误差与有效数 设2 21 1N N j S j ==-∑ ，其精确值为1311221N N ??-- ? +?? 。 （1）编制按从大到小的顺序222 111 21311 N S N = +++---，计算N S 的通用程序。 （2）编制按从小到大的顺序2221111(1)121N S N N =+++----，计算N S 的通用程序。 （3）按两种顺序分别计算210S ，410S ，610S ，并指出有效位数。（编制程序时用单精度） （4）通过本上机题，你明白了什么？ 解： （1） #include void main() { float n,i,s; printf("please input n="); scanf("%f",&n); for(i=2,s=0;i<=n;s+=1/(i*i-1),i++); printf("s=%f\n",s); } （2） #include void main() { float n,i,s; printf("please input n="); scanf("%f",&n); for(i=n,s=0;i>=2;s+=1/(i*i-1),i--); printf("s=%f\n",s); } （3）按从大到小顺序：210S =0.740049 有效位数6位 410S =0.749852 有效位数3位 610S =0.749852 有效位数3位 按从小到大顺序： 210S =0.740050 有效位数5位

410S =0.749900 有效位数6位 610S =0.749999 有效位数6位 (4)通过上述实验数据可以看出此次算法使用从小到大的顺序进行得到的数据相对而言更精确，可以得到这样的启示：在计算数值时，要先分析不同算法对结果的影响，避免大数吃小数的现象，找出能得到更精确的结果的算法。 第二章 （上机题）Newton 迭代法 （1）给定初值0x 及容许误差ε，编制Newton 法解方程()0f x =根的通用程序。 （2）给定方程3 ()/30f x x x =-=，易知其有三个根1x *=，20x * =，3x * = 1．由Newton 方法的局部收敛性可知存在0δ>，当0(,)x δδ∈-时，Newton 迭代序列收敛于根2x * 。试确定尽可能大的δ。 2．试取若干初始值，观察当0(,1)x ∈-∞-，(1,)δ--，(,)δδ-，(,1)δ，(1,)∞时Newton 序列是否收敛以及收敛于哪一个根。 （3）通过本上机题，你明白了什么？ 解： （1）#include #include #define eps 0.000001 float f(float x) { float f; f=x*x*x/3-x; return(f); } float df(float x) { float df; df=x*x-1; return (df); } void main(void) { float x0,x1,a;