计算一长一短两个序列的线性卷积时，如果直接使用 FFT 以及 IFFT 进行快速变换，则需要在较短的序列后面补很多个零，增加了计算量。因此在此种情况下，一般采用重叠相加法和重叠保留法来计算，可以有效减少计算量，具体原理此处不涉及，只提供这两种方法的matlab程序，供参考。
overlap_add 重叠相加法的程序
overlap_save 重叠保留法的程序
lin_conv 用 FFT 和 IFFT 快速计算线性卷积的程序，被以上两个程序调用
建议：将这三段程序保存成分开的三个 .m 文件
%% 重叠相加法，重叠保留法实现长短序列的线性卷积
function y=overlap_add(x,h,M,choice)
% x为较长的输入序列，h为较短的系统函数，M为分段大小
% choice为'add'或'save',选择重叠相加法或重叠保留法
%---------前期处理-------------------
N=length(h); % 序列h(n)的长度
Lx=length(x); % 序列x(n)的长度
if N>M % 算法要求 N<=M
M=N+1;
end
L=M+N-1; % 用每段循环卷积计算线性卷积所需
点数
T=ceil(Lx/M); % 分段数，ceil向上取整
t=zeros(1,N-1); % 缓存序列初始化
x=[x,zeros(1,(T+1)*M-Lx)]; % 最后一个不足M的分段补零
y=zeros(1,(T+1)*M); % 生成输出序列y(n)
%-------------核心算法--------------------
for i=0:T
xi=i*M+1;
x_seg=x(xi:xi+M-1); % 低点数计算时的分段x(n)
y_seg=lin_conv(x_seg,h,L); % 计算 i 分段和 h 的循环线性卷积
y_seg(1:N-1)=y_seg(1:N-1)+t(1:N-1); % 与前一段卷积的后N-1位重叠相
加
t(1:N-1)=y_seg(M+1:L); % 缓存序列更新
y(xi:xi+M-1)=y_seg(1:M); % 每卷积一段输出M个点
end
y=y(1:Lx+N-1) ; % 最终输出序列
end
%% 重叠保留法实现长短序列的线性卷积
function y=overlap_save(x,h,M)
% x为较长的输入序列，h为较短的系统函数，M为分段大小
% choice为'add'或'save',选择重叠相加法或重叠保留法
%---------前期处理-------------------
N=length(h); % 序列h(n)的长度
Lx=length(x); % 序列x(n)的长度
if N>M % 算法要求 N<=M
M=N+1;
end
L=M+N-1; % 用每段循环卷积计算线性卷积所需点数
T=ceil(Lx/M); % 分段数，ceil向上取整
t=zeros(1,N-1); % 缓存序列初始化
x=[x,zeros(1,(T+1)*M-Lx)]; % 最后一个不足M的分段补零
y=zeros(1,(T+1)*M); % 生成输出序列y(n)
%-------------核心算法--------------------
for i=0:T;
xi=i*M+1;
x_seg=[t,x(xi:xi+M-1)]; % 低点数计算时的分段x(n)
t=x_seg(M+1:M+N-1); % 缓存序列更新
y_seg=lin_conv(x_seg,h,L); % 线性卷积
y(xi:xi+M-1)=y_seg(N:M+N-1); % 每卷积一段输出M个点
end
y=y(1:Lx+N-1) ; % 最终输出序列
end
%% 利用FFT循环卷积计算线性卷积的函数function y=lin_conv(x1,x2,L)
% 利用FFT循环卷积计算线性卷积
y1=fft(x1,L);
y2=fft(x2,L);
Yk=y1.*y2;
y=ifft(Yk);
end
实验十七牛顿迭代法 【实验目的】 1.了解牛顿迭代法的基本概念。 2.了解牛顿迭代法的收敛性和收敛速度。 3.学习、掌握MATLAB软件的有关命令。 【实验内容】 用牛顿迭代法求方程的近似根，误差不超过。 3210 ++-=3 10- x x x 【实验准备】 1.牛顿迭代法原理 2.牛顿迭代法的几何解析 3.牛顿迭代法的收敛性 4.牛顿迭代法的收敛速度 5.迭代过程的加速 6.迭代的MATLAB命令 MATLAB中主要用for,while等控制流命令实现迭代。 【实验重点】 1.牛顿迭代法的算法实现 2.牛顿迭代法收敛性和收敛速度 【实验难点】 1.牛顿迭代法收敛性和收敛速度 【实验方法与步骤】 练习1用牛顿迭代法求方程在x=0.5附近的近似 3210 ++-= x x x
根，误差不超过。 310-牛顿迭代法的迭代函数为 322()1()()321 f x x x x g x x x f x x x ++-=-=-'++相应的MATLAB 代码为 >>clear; >>x=0.5; >>for i=1:3 >>x=x-(x^3+x^2+x-1)/(3*x^2+2*x+1) >>end 可算的迭代数列的前3项0.5455，0.5437，0.5437。经三次迭代，就大大超过了精度要求。 练习2 用牛顿迭代法求方程的近似正实根，由此建2(0)x a a =>立一种求平方根的计算方法。 由计算可知，迭代格式为，在实验12的练习4中1()()2a g x x x =+已经进行了讨论。 【练习与思考】 1.用牛顿迭代法求方程的近似根。 ln 1x x =2.为求出方程的根，在区间[1,2]内使用迭代函数进行310x x --=迭代，纪录迭代数据，问迭代是否收敛？对迭代进行加速，对比加速前的数据，比较加速效果。 3.使用在不动点的泰勒公式，证明牛顿迭代法收敛原理。*x
牛顿插值法 插值法是利用函数f (x)在某区间中若干点的函数值，作出适当的特定函数，在这些点上取已知值，在区间的其他点上用这特定函数的值作为函数f (x)的近似值。如果这特定函数是多项式，就称它为插值多项式。当插值节点增减时全部插值基函数均要随之变化，这在实际计算中很不方便。为了克服这一缺点，提出了牛顿插值。牛顿插值通过求各阶差商，递推得到的一个公式： f(x)=f[x0]+f[x0,x1](x-x0)+f[x0,x1,x2](x-x0)(x-x1)+...f[x0,...xn](x-x0 )...(x-xn-1)+Rn(x)。 插值函数 插值函数的概念及相关性质[1] 定义：设连续函数y-f(x) 在区间[a,b]上有定义，已知在n+1个互异的点 x0,x1,…xn上取值分别为y0,y1,…yn （设a≤ x1≤x2……≤xn≤b)。若在函数类中存在以简单函数P(x) ，使得P(xi)=yi,则称P(x) 为f(x)的插值函数. 称x1,x2,…xn 为插值节点，称[a,b]为插值区间。 定理：n次代数插值问题的解存在且唯一。
牛顿插值法C程序 程序框图#include Several Typical Interpolation in Matlab Lagrange Interpolation Supposing： If x=175, while y=? Solution： Lagrange Interpolation in Matlab: function y=lagrange(x0,y0,x); n=length(x0);m=length(x); for i=1:m z=x(i); s=0.0; for k=1:n p=1.0; for j=1:n if j~=k p=p*(z-x0(j))/(x0(k)-x0(j)); end end s=p*y0(k)+s; end y(i)=s; end input： x0=[144 169 225] y0=[12 13 15] y=lagrange(x0,y0,175) obtain the answer： x0 = 144 169 225 y0 = 12 13 15 y = 13.2302 Spline Interpolation Solution ： Input x=[ 1 4 9 6]；y=[ 1 4 9 6]；x=[ 1 4 9 6]；pp=spline(x,y) pp = form: 'pp' breaks: [1 4 6 9] coefs: [3x4 double] pieces: 3 order: 4 dim: 1 output : pp.coefs ans = -0.0500 0.5333 -0.8167 1.0000 -0.0500 0.0833 1.0333 2.0000 -0.0500 -0.2167 0.7667 4.0000 It shows the coefficients of cubic spline polynomial , so: S （x ）=, 169,3)9(1484.0)9(0063.0)9(0008.0,94,2)4(2714.0)4(0183.0)4(0008 .0, 41,1)1(4024.0)1(0254.0)1(0008.0232 3 23≥≤+-+---≥≤+-+---≥≤+-+---x x x x x x x x x x x x Newton’s Interpolation Resolve 65 Solution: Newton’s Interpolation in matlab ： function yi=newint(x,y,xi); n=length(x); ny=length(y); if n~=ny error end Y=zeros(n);Y(:,1)=y'; 解线性方程组的迭代法 1.rs里查森迭代法求线性方程组Ax=b的解 function[x,n]=rs(A,b,x0,eps,M) if(nargin==3) eps=1.0e-6;%eps表示迭代精度 M=10000;%M表示迭代步数的限制值elseif(nargin==4) M=10000; end I=eye(size(A)); n=0; x=x0; tol=1; %迭代过程 while(tol>eps) x=(I-A)*x0+b; n=n+1;%n为最终求出解时的迭代步数tol=norm(x-x0); x0=x; if(n>=M) disp('Warning:迭代次数太多，可能不收敛！'); return; end end 2.crs里查森参数迭代法求线性方程组Ax=b的解 function[x,n]=crs(A,b,x0,w,eps,M) if(nargin==4) eps=1.0e-6;%eps表示迭代精度 M=10000;%M表示迭代步数的限制值 elseif(nargin==5) M=10000; end I=eye(size(A)); n=0; x=x0; tol=1; %迭代过程 while(tol>eps) x=(I-w*A)*x0+w*b; n=n+1;%n为最终求出解时的迭代步数tol=norm(x-x0); x0=x; if(n>=M) disp('Warning:迭代次数太多，可能不收敛！'); return; end end 3.grs里查森迭代法求线性方程组Ax=b的解 function[x,n]=grs(A,b,x0,W,eps,M) if(nargin==4) eps=1.0e-6;%eps表示迭代精度 M=10000;%M表示迭代步数的限制值 elseif(nargin==5) M=10000; end I=eye(size(A)); n=0; x=x0; tol=1;%前后两次迭代结果误差 %迭代过程 while(tol>eps) x=(I-W*A)*x0+W*b;%迭代公式 n=n+1;%n为最终求出解时的迭代步数tol=norm(x-x0); x0=x; if(n>=M) disp('Warning:迭代次数太多，可能不收敛！'); return; end end 4.jacobi雅可比迭代法求线性方程组Ax=b的解 function[x,n]=jacobi(A,b,x0,eps,varargin) if nargin==3 eps=1.0e-6; M=200; elseif nargin<3 error return elseif nargin==5 M=varargin{1}; end D=diag(diag(A));%求A的对角矩阵 L=-tril(A,-1);%求A的下三角阵 Matlab实现数值分析插值及积分 摘要： 数值分析(numerical analysis)是研究分析用计算机求解数学计算问题的数值计算方法及其理论的学科，是数学的一个分支，它以数字计算机求解数学问题的理论和方法为研究对象。在实际生产实践中，常常将实际问题转化为数学模型来解决，这个过程就是数学建模。学习数值分析这门课程可以让我们学到很多的数学建模方法。 分别运用matlab数学软件编程来解决插值问题和数值积分问题。题目中的要求是计算差值和积分，对于问题一，可以分别利用朗格朗日插值公式，牛顿插值公式，埃特金逐次线性插值公式来进行编程求解，具体matlab代码见正文。编程求解出来的结果为：=+。 其中Aitken插值计算的结果图如下： 对于问题二，可以分别利用复化梯形公式，复化的辛卜生公式，复化的柯特斯公式编写程序来进行求解，具体matlab代码见正文。编程求解出来的结果为： 0.6932 其中复化梯形公式计算的结果图如下： 问题重述 问题一：已知列表函数 表格 1 分别用拉格朗日，牛顿，埃特金插值方法计算。 问题二：用复化的梯形公式，复化的辛卜生公式，复化的柯特斯公式计算积分，使精度小于5。 问题解决 问题一：插值方法 对于问题一，用三种差值方法：拉格朗日，牛顿，埃特金差值方法来解决。 一、拉格朗日插值法： 拉格朗日插值多项式如下： 首先构造1+n 个插值节点n x x x ,,,10 上的n 插值基函数，对任一点i x 所对应的插值基函数 )(x l i ，由于在所有),,1,1,,1,0(n i i j x j +-=取零值，因此)(x l i 有因子 )())(()(110n i i x x x x x x x x ----+- 。又因)(x l i 是一个次数不超过n 的多项式，所以只 可能相差一个常数因子，固)(x l i 可表示成： )())(()()(110n i i i x x x x x x x x A x l ----=+- 利用1)(=i i x l 得： Matlab编程实现主成分分析 .程序结构及函数作用 在软件Matlab中实现主成分分析可以采取两种方式实现：一是通过编程来实现；二是直接调用Matlab种自带程序实现。下面主要主要介绍利用Matlab的矩阵计算功能编程实现主成分分析。 1程序结构 2函数作用 Cwstd.m——用总和标准化法标准化矩阵 Cwfac.m——计算相关系数矩阵；计算特征值和特征向量；对主成分进行排序；计算各特征值贡献率；挑选主成分（累计贡献率大于85%），输出主成分个数；计算主成分载荷 Cwscore.m——计算各主成分得分、综合得分并排序 Cwprint.m——读入数据文件；调用以上三个函数并输出结果 3.源程序 3.1 cwstd.m总和标准化法标准化矩阵 %cwstd.m,用总和标准化法标准化矩阵 function std=cwstd(vector) cwsum=sum(vector,1); %对列求和 [a,b]=size(vector); %矩阵大小,a为行数,b为列数 for i=1:a for j=1:b std(i,j)= vector(i,j)/cwsum(j); end end 3.2 cwfac.m计算相关系数矩阵 %cwfac.m function result=cwfac(vector); fprintf('相关系数矩阵:\n') std=CORRCOEF(vector) %计算相关系数矩阵 fprintf('特征向量(vec)及特征值(val)：\n') [vec,val]=eig(std) %求特征值(val)及特征向量(vec) newval=diag(val) ; [y,i]=sort(newval) ; %对特征根进行排序，y为排序结果，i为索引fprintf('特征根排序：\n') for z=1:length(y) newy(z)=y(length(y)+1-z); end fprintf('%g\n',newy) rate=y/sum(y); fprintf('\n贡献率：\n') newrate=newy/sum(newy) sumrate=0; newi=[]; for k=length(y):-1:1 sumrate=sumrate+rate(k); newi(length(y)+1-k)=i(k); if sumrate>0.85 break; end end %记下累积贡献率大85%的特征值的序号放入newi中fprintf('主成分数：%g\n\n',length(newi)); fprintf('主成分载荷：\n') for p=1:length(newi) 拉格朗日插值的调用函数 function y=lagrange(x0,y0,x) n=length(x0);m=length(x); for i=1:m z=x(i); L=0.0; for j=1:n T=1.0; for k=1:n if k~=j T=T*(z-x0(k))/(x0(j)-x0(k)); end end L=T*y0(j)+L; end y(i)=L; end 四个图在一起： x=[-1:0.05:1]; y=1./(1+25*x.^2); x0=[-1:0.4:1]; y0=1./(1+25*x0.^2); y1=lagrange(x0,y0,x); x0=[-1:0.2:1]; y0=1./(1+25*x0.^2); y2=lagrange(x0,y0,x); x0=[-1:0.1:1]; y0=1./(1+25*x0.^2); y3= lagrange(x0,y0,x); plot(x,y,'-r') hold on plot(x,y1,'-b',x,y2,'-r',x,y3,'-r') l5和fx在一起： x=[-1:0.05:1]; y=1./(1+25*x.^2); x0=[-1:0.4:1]; y0=1./(1+25*x0.^2); y1=lagrange(x0,y0,x); plot(x,y,'-r') hold on plot(x,y1,'-b') l10和fx在一起： x=[-1:0.05:1]; y=1./(1+25*x.^2); x0=[-1:0.2:1]; y0=1./(1+25*x0.^2); y2= lagrange(x0,y0,x); plot(x,y,'-r') hold on plot(x,y2,'-b') l20和fx在一起： x=[-1:0.05:1]; y=1./(1+25*x.^2); x0=[-1:0.1:1]; y0=1./(1+25*x0.^2); y3= lagrange(x0,y0,x); plot(x,y,'-r') hold on plot(x,y3,'-b') 层次分析法建模 层次分析法（AHP－Analytic Hierachy process）---- 多目标决策方法 70 年代由美国运筹学家T·L·Satty提出的，是一种定性与定量分析相结合的多目标决策分析方法论。吸收利用行为科学的特点，是将决策者的经验判断给予量化，对目标（因素）结构复杂而且缺乏必要的数据情况下，採用此方法较为实用，是一种系统科学中，常用的一种系统分析方法，因而成为系统分析的数学工具之一。 传统的常用的研究自然科学和社会科学的方法有： 机理分析方法：利用经典的数学工具分析观察的因果关系； 统计分析方法：利用大量观测数据寻求统计规律，用随机数学方法描述（自然现象、 社会现象）现象的规律。 基本内容：（1）多目标决策问题举例AHP建模方法 （2）AHP建模方法基本步骤 （3）AHP建模方法基本算法 （3）AHP建模方法理论算法应用的若干问题。 参考书：1、姜启源，数学模型（第二版，第9章；第三版，第8章），高等教育出版社 2、程理民等，运筹学模型与方法教程，（第10章），清华大学出版社 3、《运筹学》编写组，运筹学（修订版），第11章，第7节，清华大学出版社 一、问题举例： A．大学毕业生就业选择问题 获得大学毕业学位的毕业生，“双向选择”时，用人单位与毕业生都有各自的选择标准和要求。就毕业生来说选择单位的标准和要求是多方面的，例如： ①能发挥自己的才干为国家作出较好贡献（即工作岗位适合发挥专长）； ②工作收入较好（待遇好）； ③生活环境好（大城市、气候等工作条件等）； ④单位名声好（声誉-Reputation）； ⑤工作环境好（人际关系和谐等） ⑥发展晋升（promote, promotion）机会多（如新单位或单位发展有后劲）等。 问题：现在有多个用人单位可供他选择，因此，他面临多种选择和决策，问题是他将如何作出决策和选择？——或者说他将用什么方法将可供选择的工作单位排序？ 一、实验目的及题目 1.1 实验目的： （1）学会用高斯列主元消去法，LU 分解法，Jacobi 迭代法和Gauss-Seidel 迭代法解线性方程组。 （2）学会用Matlab 编写各种方法求解线性方程组的程序。 1.2 实验题目： 1. 用列主元消去法解方程组： 1241234 123412343421233234x x x x x x x x x x x x x x x ++=??+-+=??--+=-??-++-=? 2. 用LU 分解法解方程组,Ax b =其中 4824012242412120620266216A --?? ?- ?= ? ?-??，4422b ?? ? ?= ?- ?-?? 3. 分别用Jacobi 迭代法和Gauss-Seidel 迭代法求解方程组： 123234 1231234102118311210631125x x x x x x x x x x x x x -+=-??-+=-??-+=??-+-+ =? 二、实验原理、程序框图、程序代码等 2.1实验原理 2.1.1高斯列主元消去法的原理 Gauss 消去法的基本思想是一次用前面的方程消去后面的未知数，从而将方程组化为等价形式： 1111221122222n n n n nn n n b x b x b x g b x b x g b x g +++=??++=????= ? 这个过程就是消元，然后再回代就好了。具体过程如下： 对于1,2, ,1k n =-，若() 0,k kk a ≠依次计算 ()() (1)()()(1)()()/,,1, ,k k ik ik kk k k k ij ij ik kj k k k i i ik k m a a a a m a b b m b i j k n ++==-=-=+ 然后将其回代得到： ()() ()()()1/()/,1,2,,1 n n n n nn n k k k k k kj j kk j k x b a x b a x a k n n =+?=??=-=--? ? ∑ 以上是高斯消去。 但是高斯消去法在消元的过程中有可能会出现() 0k kk a =的情况，这时消元就无法进行了，即使主元数() 0,k kk a ≠但是很小时，其做除数，也会导致其他元素数量级的严重增长和舍入误差的扩散。因此，为了减少误差，每次消元选取系数矩阵的某列中绝对值最大的元素作为主元素。然后换行使之变到主元位置上，再进行销元计算。即高斯列主元消去法。 2.1.2直接三角分解法（LU 分解）的原理 先将矩阵A 直接分解为A LU =则求解方程组的问题就等价于求解两个三角形方程组。 直接利用矩阵乘法，得到矩阵的三角分解计算公式为： 1111111 11 1,1,2,,/,2,,,,,1,,,2,3, ()/,1,2, ,i i i i k kj kj km mj m k ik ik im mk kk m u a i n l a u i n u a l u j k k n k n l a l u u i k k n k n -=-===?? ==?? =-=+??=??=-=++≠?? ∑∑且 由上面的式子得到矩阵A 的LU 分解后，求解Ux=y 的计算公式为 11 111,2,3,/()/,1,2, ,1 i i i ij j j n n nn n i i ij j ii j i y b y b l y i n x y u x y u x u i n n -==+=??? =-=?? =??? =-=--?? ∑∑ 以上为LU 分解法。 数值分析作业 姓名王建忠 学号132080202006 学院能源与动力工程 专业机械电子工程 2013年12月16日 https://www.doczj.com/doc/dd16431413.html,grange插值多项式程序 function f=nalagr(x,y,xx) %x为节点向量 %y为节点函数值 %xx是插值点 syms s if(length(x)==length(y)) n=length(x); else disp('x和y的维数不相等!'); return; end f=0.0; for(i=1:n) l=y(i); for(j=1:i-1) l=l*(s-x(j))/(x(i)-x(j)); end; for(j=i+1:n) l=l*(s-x(j))/(x(i)-x(j));%计算拉格朗日基函数end; f=f+l;%计算拉格朗日插值函数 simplify(f); if(i==n) if(nargin==3) f=subs(f,'s');%计算插值点的函数值else f=collect(f);%将插值多项式展开 f=vpa(f,6);%将插值多项式的系数化成6位精度的小数 end end end >>x=[-2,-1,0,1];%已知节点向量y=[3,1,1,6];%节点函数值向量 f=nalagr(x,y) f= 0.5*s^3+ 2.5*s^2+ 2.0*s+ 1.0 >>f=nalagr(x,y,0) f=1 >> 2.牛顿插值多项式程序 function[p2,z]=newTon(x,y,t) %输入参数中x,y为元素个数相等的向量，t为待估计的点，可以为数字或向量。%输出参数中p2为所求得的牛顿插值多项式，z为利用多项式所得的t的函数值。 n=length(x); chaS(1)=y(1); for i=2:n x1=x;y1=y; x1(i+1:n)=[]; y1(i+1:n)=[]; n1=length(x1); s1=0; for j=1:n1 t1=1; for k=1:n1 if k==j continue; else t1=t1*(x1(j)-x1(k)); end end s1=s1+y1(j)/t1; end chaS(i)=s1; end b(1,:)=[zeros(1,n-1)chaS(1)]; cl=cell(1,n-1); for i=2:n u1=1; for j=1:i-1 u1=conv(u1,[1-x(j)]); cl{i-1}=u1; end cl{i-1}=chaS(i)*cl{i-1}; b(i,:)=[zeros(1,n-i),cl{i-1}]; end p2=b(1,:); for j=2:n p2=p2+b(j,:); end if length(t)==1 rm=0; 要求： 下面分别使用雅克比迭代法和高斯-赛德尔迭代法求一个方程组的近似解用的线性方程组是按实验要求给的： 7*x1+x2+2*x3=10 x1+8*x2+2*x3=8 2*x1+2*x2+9*x3=6 雅克比迭代法的matlab代码：（老师写的） A=[7,1,2;1,8,2;2,2,9]; b=[10;8;6]; if(any(diag(A))==0) error('error,pause') end eps=input('误差限eps='); N=input('迭代次数N='); D=diag(diag(A)); B=inv(D)*(D-A); f=inv(D)*b; K=0; x0=zeros(size(b)); while 1 x1=B*x0+f K=K+1; fprintf('第-次迭代的近似解为',K) disp(x1'); if norm(x1-x0,inf) A=[7,1,2;1,8,2;2,2,9]; b=[10;8;6]; if(all(diag(A))==0) error('error,pause') end eps=input('误差限eps='); N=input('迭代次数N='); D=diag(diag(A)); B=inv(D)*(D-A); f=inv(D)*b; K=0; x0=zeros(size(b)); x00=x0; while 1 x11=B*x0+f; x00(1,1)=x11(1,1); x12=B*x00+f; x00(2,1)=x12(2,1); x13=B*x00+f; x00(3,1)=x13(3,1); x1=x00 K=K+1; fprintf('第-次迭代的近似解为',K) disp(x1'); if norm(x1-x0,inf) Matlab中插值函数汇总和使用说明 MATLAB中的插值函数为interp1，其调用格式 为： yi= interp1(x,y,xi,'method') 其中x，y为插值点，yi为在被插值点xi处的插值结果；x,y为向量， 'method'表示采用的插值方法，MATLAB提供的插值方法有几种： 'method'是最邻近插值， 'linear'线性插值； 'spline'三次样条插值； 'cubic'立方插值．缺省时表示线性插值 注意：所有的插值方法都要求x是单调的，并且xi不能够超过x的范围。 例如：在一天24小时内，从零点开始每间隔2小时测得的环境温度数据分别为 12，9，9，10，18 ，24，28，27，25，20，18，15，13， 推测中午12点（即13点）时的温度． x=0:2:24; y=[12 9 9 10 18 24 28 27 25 20 18 15 13]; a=13; y1=interp1(x,y,a,'spline') 结果为： 27.8725 若要得到一天24小时的温度曲线，则： xi=0:1/3600:24; yi=interp1(x,y,xi, 'spline'); plot(x,y,'o' ,xi,yi) 命令1 interp1 功能一维数据插值（表格查找）。该命令对数据点之间计算内插值。它找出一元函数f(x)在中间点的数值。其中函数f(x)由所给数据决定。 x：原始数据点 Y：原始数据点 xi：插值点 Yi：插值点 格式 (1)yi = interp1(x,Y,xi) 返回插值向量yi，每一元素对应于参量xi，同时由向量x 与Y 的内插值决定。参量x 指定数据Y 的点。 若Y 为一矩阵，则按Y 的每列计算。yi是阶数为length(xi)*size(Y,2)的输出矩阵。(2)yi = interp1(Y,xi) 假定x=1:N，其中N 为向量Y 的长度，或者为矩阵Y 的行数。 (3)yi = interp1(x,Y,xi,method) 用指定的算法计算插值： ’nearest’：最近邻点插值，直接完成计算； ’linear’：线性插值（缺省方式），直接完成计算； ’spline’：三次样条函数插值。对于该方法，命令interp1 调用函数spline、ppval、mkpp、umkpp。这些命令生成一系列用于分段多项式操作的函数。命令spline 用它们执行三次样条函数插值； ’pchip’：分段三次Hermite 插值。对于该方法，命令interp1 调用函数pchip，用于对 1 1.微分方程的定义 对于duffing 方程03 2 =++x x x ω ，先将方程写作??? --==3 1122 21x x x x x ω function dy=duffing(t,x) omega=1;%定义参数 f1=x(2); f2=-omega^2*x(1)-x(1)^3; dy=[f1;f2]; 2.微分方程的求解 function solve (tstop) tstop=500;%定义时间长度 y0=[0.01;0];%定义初始条件 [t,y]=ode45('duffing',tstop,y0,[]); function solve (tstop) step=0.01;%定义步长 y0=rand(1,2);%随机初始条件 tspan=[0:step:500];%定义时间范围 [t,y]=ode45('duffing',tspan,y0); 3.时间历程的绘制 时间历程横轴为t ，纵轴为y ，绘制时只取稳态部分。 plot(t,y(:,1));%绘制y 的时间历程 xlabel('t')%横轴为t ylabel('y')%纵轴为y grid;%显示网格线 2 axis([460 500 -Inf Inf])%图形显示范围设置 4.相图的绘制 相图的横轴为y ，纵轴为dy/dt ，绘制时也只取稳态部分。红色部分表示只取最后1000个点。 plot(y(end-1000:end ,1),y(end-1000:end ,2));%绘制y 的时间历程 xlabel('y')%横轴为y ylabel('dy/dt')%纵轴为dy/dt grid;%显示网格线 5.Poincare 映射的绘制 对于不同的系统，Poincare 截面的选取方法也不同 对于自治系统一般每过其对应线性系统的固有周期，截取一次 对于非自治系统，一般每过其激励的周期，截取一次 例程：duffing 方程03 2=++x x x ω 的poincare 映射 function poincare(tstop) global omega; omega=1; T=2*pi/omega;%线性系统的周期或激励的周期 step=T/100;%定义步长为T/100 y0=[0.01;0];%初始条件 tspan=[0:step:100*T];%定义时间范围 [t,y]=ode45('duffing',tspan,y0); for i=5000:100:10000%稳态过程每个周期取一个点 plot(y(i,1),y(i,2),'b.'); hold on;% 保留上一次的图形 end xlabel('y');ylabel('dy/dt'); 题目一：多项式插值 某气象观测站在8：00（AM ）开始每隔10分钟对天气作如下观测，用三次多项式插值函数（Newton ）逼近如下曲线，插值节点数据如上表，并求出9点30分该地区的温度（x=10）。 二、数学原理 假设有n+1个不同的节点及函数在节点上的值（x 0，y 0），……（x n ，y n ），插值多项式有如下形式： ）（） ）（（）（）（）（n 10n 102010n x -x )(x -x x -x x P x x x x x x -??-+??+-++=αααα （1） 其中系数i α（i=0,1,2……n ）为特定系数，可由插值样条i i n y x P =） （（i=0,1,2……n ）确定。 根据均差的定义，把x 看成[a,b]上的一点,可得 f(x)= f （0x ）+f[10x x ，]（0x -x ） f[x, 0x ]= f[10x x ，]+f[x,10x x ，] （1x -x ） …… f[x, 0x ，…x 1-n ]= f[x, 0x ，…x n ]+ f[x, 0x ，…x n ](x-x n ) 综合以上式子，把后一式代入前一式，可得到： f(x)= f[0x ]+f[10x x ，]（0x -x ）+ f[210x x x ，，]（0x -x ）（1x -x ）+ …+ f[x, 0x ，…x n ]（0x -x ）…(x-x 1-n )+ f[x, 0x ，…x n ,x ]） （x 1n +ω= N n （x ）+） （x n R 其中 N n （x ）= f[0x ]+f[10x x ，]（0x -x ）+ f[210x x x ，，]（0x -x ）（1x -x ）+ …+ f[x, 0x ，…x n ]（0x -x ）…(x-x 1-n ) （2） ）（x n R = f(x)- N n （x ）= f[x, 0x ， (x) n ,x ]） （x 1n +ω （3） ） （x 1n +ω=（0x -x ）…(x-x n ) Newton 插值的系数i α（i=0,1,2……n ）可以用差商表示。一般有 f k =α[k 10x x x ??，] （k=0,1,2，……，n ） （4） 把（4）代入（1）得到满足插值条件N ）（） （i i n x f x =（i=0,1,2，……n ）的n 次Newton 插值多项式 N n （x ）=f （0x ）+f[10x x ，]（1x -x ）+f[210x x x ，，]（1x -x ）（2x -x ）+……+f[n 10x x x ??，]（1x -x ）（2x -x ）…（1-n x -x ）. 其中插值余项为： ） （）！ （） （）（）（）（x 1n f x N -x f x R 1n 1 n n +++==ωξ ξ介于k 10x x x ??，之间。 三、程序设计 function [y,A,C,L]=newdscg(X,Y,x,M) % y 为对应x 的值，A 为差商表，C 为多项式系数，L 为多项式 % X 为给定节点，Y 为节点值，x 为待求节点 n=length(X); m=length(x); % n 为X 的长度 for t=1:m 实验一非线性方程的数值解法(一) 信息与计算科学金融崔振威201002034031一、实验目的： 熟悉二分法和简单迭代法的算法实现。 二、实验内容： 教材P40 2.1.5 三、实验要求 1根据实验内容编写二分法和简单迭代法的算法实现 2简单比较分析两种算法的误差 3试构造不同的迭代格式，分析比较其收敛性 (一)、二分法程序： function ef=bisect(fx,xa,xb ,n, delta) % fx是由方程转化的关于x的函数，有fx=0。 % xa解区间上限 % xb解区间下限 % n最多循环步数，防止死循环。 %delta为允许误差 x=xa;fa=eval(fx); x=xb;fb=eval(fx); disp(' [ n xa xb xc fc ]'); for i=1: n xc=(xa+xb)/2;x=xc;fc=eval(fx); X=[i,xa,xb,xc,fc]; disp(X), if fc*fa<0 xb=xc; else xa=xc; end if (xb-xa) k=0; while abs(x-xO)>eps & k MATLAB程序设计期中作业 ——编程实现牛顿插值 成员:刘川（P091712797）签名＿＿＿＿＿ 汤意（P091712817）签名＿＿＿＿＿ 王功贺（P091712799）签名＿＿＿＿＿ 班级:2009信息与计算科学 学院:数学与计算机科学学院 日期:2012年05月02日 牛顿插值的算法描述及程序实现 一：问题说明 在我们的实际应用中，通常需要解决这样的问题，通过一些已知的点及其对应的值，去估算另外一些点的值，这些数据之间近似服从一定的规律，于是，这就引入了插值法的思想。 插值法是利用函数f (x)在某区间中若干点的函数值，作出适当的特定函数，在这些点上取已知值，在区间的其他点上用这特定函数的值作为函数f (x)的近似值。如果这特定函数是多项式，就称它为插值多项式。利用插值基函数很容易得到拉格朗日插值多项式，公式结构紧凑，在理论分析中甚为方便，但当插值节点增减时全部插值基函数均要随之变化，整个公式也将发生变化，这在实际计算中是很不方便的，为了克服这一缺点，提出了牛顿插值。 二：算法分析 newton 插值多项式的表达式如下： 010011()()()()()n n n N x c c x x c x x x x x x -=+-+???+--???- 其中每一项的系数c i 的表达式如下： 12011010 [,,,][,,,] [,,,]i i i i i f x x x f x x x c f x x x x x -???-???=???= - 即为f (x)在点01,,,i x x x ???处的i 阶差商，（[]()i i f x f x =，1,2,,i n = ），由差商01[,,,]i f x x x ???的性质可知： () 010 1 [,,,]()i i i j j k j k k j f x x x f x x x ==≠???=-∑∏ 牛顿插值的程序实现方法： 第一步：计算[][][][]001012012,,,,,,,n f x f x x f x x x f x x x x 、、、 、。 第二步：计算牛顿插值多项式中01[,,,]i f x x x ???011()()()i x x x x x x ---???-，1,2,,i n = ，得到n 个多项式。 6.1 插值问题及其误差 6.1.2 与插值有关的MATLAB 函数 (一) POLY2SYM函数 调用格式一：poly2sym (C) 调用格式二：f1=poly2sym(C,'V') 或f2=poly2sym(C, sym ('V') ), (二) POLYVAL函数 调用格式：Y = polyval(P,X) (三) POLY函数 调用格式：Y = poly (V) (四) CONV函数 调用格式：C =conv (A, B) 例 6.1.2求三个一次多项式、和的积.它们的零点分别依次为0.4,0.8,1.2. 解我们可以用两种MATLAB程序求之. 方法1如输入MATLAB程序 >> X1=[0.4,0.8,1.2]; l1=poly(X1), L1=poly2sym (l1) 运行后输出结果为 l1 = 1.0000 - 2.4000 1.7600 -0.3840 L1 = x^3-12/5*x^2+44/25*x-48/125 方法2如输入MATLAB程序 >> P1=poly(0.4);P2=poly(0.8);P3=poly(1.2); C =conv (conv (P1, P2), P3) , L1=poly2sym (C) 运行后输出的结果与方法1相同. (五) DECONV 函数 调用格式：[Q,R] =deconv (B,A) (六) roots(poly(1:n))命令 调用格式：roots(poly(1:n)) (七) det(a*eye(size (A)) - A)命令 调用格式：b=det(a*ey e(size (A)) - A) 6.2 拉格朗日（Lagrange）插值及其MATLAB程序 6.2.1 线性插值及其MATLAB程序 例 6.2.1 已知函数在上具有二阶连续导数，，且满足条件 .求线性插值多项式和函数值，并估计其误差. 解输入程序 >> X=[1,3];Y=[1,2]; l01= poly(X(2))/( X(1)- X(2)), l11= poly(X(1))/( X(2)- X(1)), l0=poly2sym (l01),l1=poly2sym (l11), P = l01* Y(1)+ l11* Y(2), L=poly2sym (P),x=1.5; Y = polyval(P,x) 运行后输出基函数l0和l1及其插值多项式的系数向量P（略）、插值多项式L和插值Y为l0 = l1 = L = Y = -1/2*x+3/2 1/2*x-1/2 1/2*x+1/2 1.2500 输入程序 >> M=5;R1=M*abs((x-X(1))* (x-X(2)))/2matlab插值法实例
MATLAB代码 解线性方程组的迭代法
matlab实现数值分析报告插值及积分
主成分分析报告matlab程序
拉格朗日插值matlab程序
层次分析报告法及matlab程序
lu分解法、列主元高斯法、jacobi迭代法、gaussseidel法的原理及matlab程序
matlab插值法,迭代法程序
MATLAB样例之雅克比迭代法
matlab插值(详细 全面)
matlab动力学分析程序详解
matlab牛顿插值法例题与程序
二分法、简单迭代法的matlab代码实现
牛顿插值MATLAB算法
函数的插值方法及matlab程序
相关主题
文本预览