数值分析算法在matlab中的实现
插值MATLAB实现(牛顿差商插值误差龙格现象切比雪夫插值)
![插值MATLAB实现(牛顿差商插值误差龙格现象切比雪夫插值)](https://img.taocdn.com/s3/m/34a93e4d7dd184254b35eefdc8d376eeafaa1750.png)
插值MATLAB实现(牛顿差商插值误差龙格现象切比雪夫插值)插值是数值分析中的一种方法,通过已知数据点的函数值来估计函数在其他点的值。
MATLAB提供了多种方法来实现插值,包括牛顿差商插值、插值误差分析、龙格现象和切比雪夫插值。
下面将详细介绍这些方法的实现原理和MATLAB代码示例。
1.牛顿差商插值:牛顿差商插值是一种基于多项式插值的方法,其中差商是一个连续性的差分商。
该方法的优势在于可以快速计算多项式的系数。
以下是MATLAB代码示例:```matlabfunction [coeff] = newton_interpolation(x, y)n = length(x);F = zeros(n, n);F(:,1)=y';for j = 2:nfor i = j:nF(i,j)=(F(i,j-1)-F(i-1,j-1))/(x(i)-x(i-j+1));endendcoeff = F(n, :);end```该代码中,输入参数x和y分别表示已知数据点的x坐标和y坐标,返回值coeff表示插值多项式的系数。
2.插值误差分析:插值误差是指插值函数与原始函数之间的差异。
一般来说,通过增加插值节点的数量或使用更高次的插值多项式可以减小插值误差。
以下是MATLAB代码示例:```matlabfunction [error] = interpolation_error(x, y, x_eval)n = length(x);p = polyfit(x, y, n-1);y_eval = polyval(p, x_eval);f_eval = sin(pi*x_eval);error = abs(f_eval - y_eval);end```该代码中,输入参数x和y分别表示已知数据点的x坐标和y坐标,x_eval表示插值节点的x坐标,error表示插值误差。
3.龙格现象:龙格现象是插值多项式在等距插值节点上错误增长的现象。
第十章MATLAB的数值分析
![第十章MATLAB的数值分析](https://img.taocdn.com/s3/m/b1133bf0941ea76e58fa04f7.png)
• 第一个问题可归结为“已知函数在x0,x1,
– …,xn处的值,求函数在区间[x0,xn]内其它点处的值”,这 种问题适宜用插值方法解决。 – 插值问题可描述为:已知函数在x0,x1,…,xn处的值 y0,y1,…,yn,求函数p(x),使p(xi) = yi。
• 但对第二个问题不宜用插值方法,因为600米已超出所 给数据范围,用插值函数外推插值区间外的数据会 产生较大的误差。
– Q1=prctile(w,25); – Q3=prctile(w,75); – prctile( )函数实现计算样本的百分位数功能
分布形态的测定
• 只用集中趋势和离中趋势来表示所有数据,难免不 够准确。分析总体次数的分布形态有助于识别整个 总体的数量特征。总体的分布形态可以从两个角度 考虑,一是分布的对称程度,另一个是分布的高低。 前者的测定参数称为偏度或偏斜度,后者的测定参 数称为峰度。 • 峰度是掌握分布形态的另一指标,它能描述分布的 平缓或陡峭程度。如果峰度数值等于零,说明分布 为正态;若峰度数值大于零,说明分布呈陡峭状态; 若峰度数值小于零,说明分布形态趋于平缓。
– 解决第二个问题的常用方法是,根据地面到井下 500 处的 数据求出瓦斯浓度与地面到井下距离x之间的近似函数关 系f(x), 由f(x)求井下600米处的瓦斯浓度。
• 插值函数过已知点,拟合函数不一定过已知点。通 常, 插值主要用于求函数值,而拟合的主要目的是求 函数关系。当然,某些问题既可以用插值也可以用 拟合。
插值方法-概述
• 为什么需要插值?
(1) 函数关系y=f(x)没有明确的表达式
(2) y=f(x)表达式复杂,不便于研究和使用
-20 -15
沉陷量/mm 下沉方向为"+"
数值分析中求解非线性方程的MATLAB求解程序
![数值分析中求解非线性方程的MATLAB求解程序](https://img.taocdn.com/s3/m/5ab3ec56876fb84ae45c3b3567ec102de3bddf56.png)
数值分析中求解非线性方程的MATLAB求解程序1. fzero函数:fzero函数是MATLAB中最常用的求解非线性方程的函数之一、它使用了割线法、二分法和反复均值法等多种迭代算法来求解方程。
使用fzero函数可以很方便地求解单变量非线性方程和非线性方程组。
例如,要求解方程f(x) = 0,可以使用以下语法:``````2. fsolve函数:fsolve函数是MATLAB中求解多维非线性方程组的函数。
它是基于牛顿法的迭代算法来求解方程组。
使用fsolve函数可以非常方便地求解非线性方程组。
例如,要求解方程组F(x) = 0,可以使用以下语法:``````3. root函数:root函数是MATLAB中求解非线性方程组的函数之一、它采用牛顿法或拟牛顿法来求解方程组。
使用root函数可以非常方便地求解非线性方程组。
例如,要求解方程组F(x) = 0,可以使用以下语法:``````4. vpasolve函数:vpasolve函数是MATLAB中求解符号方程的函数。
它使用符号计算的方法来求解方程,可以得到精确的解。
vpasolve函数可以求解多变量非线性方程组和含有符号参数的非线性方程。
例如,要求解方程组F(x) = 0,可以使用以下语法:```x = vpasolve(F(x) == 0, x)```vpasolve函数会返回方程组的一个精确解x。
5. fsolve和lsqnonlin结合:在MATLAB中,可以将求解非线性方程转化为求解最小二乘问题的形式。
可以使用fsolve函数或lsqnonlin函数来求解最小二乘问题。
例如,要求解方程f(x) = 0,可以将其转化为最小二乘问题g(x) = min,然后使用fsolve或lsqnonlin函数来求解。
具体使用方法可以参考MATLAB官方文档。
6. Newton-Raphson法手动实现:除了使用MATLAB中的函数来求解非线性方程,还可以手动实现Newton-Raphson法来求解。
数值分析的matlab实现
![数值分析的matlab实现](https://img.taocdn.com/s3/m/19bc4ca469dc5022aaea005f.png)
第2章牛顿插值法实现参考文献:[1]岑宝俊. 牛顿插值法在凸轮曲线修正设计中的应用[J]. 机械工程师,2009,10:54-55.求牛顿插值多项式和差商的MA TLAB 主程序:function[A,C,L,wcgs,Cw]=newpoly(X,Y)n=length(X);A=zeros(n,n);A(:,1) =Y';s=0.0;p=1.0;q=1.0;c1=1.0;for j=2:nfor i=j:nA(i,j)=(A(i,j-1)-A(i-1,j-1))/(X(i)-X(i-j+1));endb=poly(X(j-1));q1=conv(q,b);c1=c1*j;q=q1;endC=A(n,n);b=poly(X(n));q1=conv(q1,b);for k=(n-1):-1:1C=conv(C,poly(X(k)));d=length(C);C(d)=C(d)+A(k,k);endL(k,:)=poly2sym(C);Q=poly2sym(q1);syms Mwcgs=M*Q/c1;Cw=q1/c1;(1)保存名为newpoly.m 的M 文件(2)输入MA TLAB 程序>> X=[242,243,249,250];>> Y=[13.681,13.526,13.098,13.095];>> [A,C,L,wcgs,Cw]=newpoly(X,Y)输出3阶牛顿插值多项式L 及其系数向量C 差商的矩阵A ,插值余项wcgs 及其)()()1(ξ+n n f x R 的系数向量Cw 。
A =13.6810 0 0 013.5260 -0.1550 0 013.0980 -0.0713 0.0120 013.0950 -0.0030 0.0098 -0.0003C =1.0e+003 *-0.0000 0.0002 -0.0551 4.7634L =- (23*x^3)/84000 + (2981*x^2)/14000 - (7757472138947345*x)/140737488355328 + 5237382665812919/1099511627776wcgs =(M*(x^4 - 984*x^3 + 363071*x^2 - 59535444*x + 3660673500))/24Cw =1.0e+008 *0.0000 -0.0000 0.0002 -0.0248 1.5253输入MATLAB程序>> x=244;>> y=- (23*x^3)/84000 + (2981*x^2)/14000 - (7757472138947345*x)/140737488355328 + 5237382665812919/1099511627776y =13.3976输入MATLAB程序>> x=[244,245,246,247,248];>> y=- (23.*x.^3)./84000 + (2981.*x.^2)./14000 - (7757472138947345.*x)./140737488355328 + 5237382665812919./1099511627776y =13.3976 13.2943 13.2143 13.1560 13.1178第4章 高斯-勒让德积分公式实现用高斯-勒让德积分公式计算dx e I x ⎰--=112221π,取代数精度为3和5,再根据截断误差公式写出误差公式,并将计算结果与精确值进行比较。
《数值分析》课程数值积分的Matlab实现问题的教学研究
![《数值分析》课程数值积分的Matlab实现问题的教学研究](https://img.taocdn.com/s3/m/5a53e0cd8bd63186bcebbc9a.png)
q a g udk q a g ( nx n ma ) u d kf ,mi, x u x G
自适应
高, 支持无穷区间积 分 , 供 thT 1 提 Asn, ’
a 一 rno 算 法 ’dT l. l Ko rd R o’ Ma[lrae t |等 xnevlot ’ o 选 项
思维能力 , 又能够培养学生应用科学计算方法和计算机技术分析
解 决 实 际 问题 的 能力 。在 我 目 , 乎 所 有 工 科 院 校 硕 i 研 究 生 都 几
开设 了《 数值分析》 课程 。 近年来 , 人们已意识到在数值分析课程的课堂教 学和实验教
学 中 引 入科 学 汁算 软 件 的重 要 性 , t b 件 已替 代 C语 言 成 为 Ma a 软 l
其 中 fn 以是 匿名 函数 、 u可 字符 串形式 的函数 、 内嵌 函数 、 M 函数文 件 的函数句 柄等 , 函数 表达 式要遵 循 “ 组运 算 ”x n 数 ,mi, x a, mi,ma 分 别表示 两个积 分变量 的 七下 限 , l xy ny x n 都为确定 的 常数 , 】 t 是绝对误 差限 , ehd 积分方法 的选择 , o m to 是 缺省方法 是 @q a, ud 还可 以选取@q al 自定义 的积分 方法 函数文件的函数 u d或
堕 塾查
NO. 7 T ME DUCAT1 I E 0N Jul y
《 数值分析》 课程数值积分的 Malb t 实现问题 的教学研究 a
戈 慈水
摘 要 : 文 根 据 数 值 分 析 课 程 的 实践 性 特 点 , 数 值 积 分 的 教 学 中融入 M aa 本 在 tb实现 问题 的教 学 , 合 问 题 教 学 法 , l 结 以提 高 学生 对 数 值 积 分 方 法 的理 解和 应 用 能 力 , 实现 本课 程 数值 积 分 内容 的 教 学 目标 。
数值分析matlab实验报告
![数值分析matlab实验报告](https://img.taocdn.com/s3/m/e2aae840a9114431b90d6c85ec3a87c241288a50.png)
数值分析matlab实验报告数值分析 Matlab 实验报告一、实验目的数值分析是研究各种数学问题数值解法的学科,Matlab 则是一款功能强大的科学计算软件。
本次实验旨在通过使用 Matlab 解决一系列数值分析问题,加深对数值分析方法的理解和应用能力,掌握数值计算中的误差分析、数值逼近、数值积分与数值微分等基本概念和方法,并培养运用计算机解决实际数学问题的能力。
二、实验内容(一)误差分析在数值计算中,误差是不可避免的。
通过对给定函数进行计算,分析截断误差和舍入误差的影响。
例如,计算函数$f(x) =\sin(x)$在$x = 05$ 附近的值,比较不同精度下的结果差异。
(二)数值逼近1、多项式插值使用拉格朗日插值法和牛顿插值法对给定的数据点进行插值,得到拟合多项式,并分析其误差。
2、曲线拟合采用最小二乘法对给定的数据进行线性和非线性曲线拟合,如多项式曲线拟合和指数曲线拟合。
(三)数值积分1、牛顿柯特斯公式实现梯形公式、辛普森公式和柯特斯公式,计算给定函数在特定区间上的积分值,并分析误差。
2、高斯求积公式使用高斯勒让德求积公式计算积分,比较其精度与牛顿柯特斯公式的差异。
(四)数值微分利用差商公式计算函数的数值导数,分析步长对结果的影响,探讨如何选择合适的步长以提高精度。
三、实验步骤(一)误差分析1、定义函数`compute_sin_error` 来计算不同精度下的正弦函数值和误差。
```matlabfunction value, error = compute_sin_error(x, precision)true_value = sin(x);computed_value = vpa(sin(x), precision);error = abs(true_value computed_value);end```2、在主程序中调用该函数,分别设置不同的精度进行计算和分析。
(二)数值逼近1、拉格朗日插值法```matlabfunction L = lagrange_interpolation(x, y, xi)n = length(x);L = 0;for i = 1:nli = 1;for j = 1:nif j ~= ili = li (xi x(j))/(x(i) x(j));endendL = L + y(i) li;endend```2、牛顿插值法```matlabfunction N = newton_interpolation(x, y, xi)n = length(x);%计算差商表D = zeros(n, n);D(:, 1) = y';for j = 2:nfor i = j:nD(i, j) =(D(i, j 1) D(i 1, j 1))/(x(i) x(i j + 1));endend%计算插值结果N = D(1, 1);term = 1;for i = 2:nterm = term (xi x(i 1));N = N + D(i, i) term;endend```3、曲线拟合```matlab%线性最小二乘拟合p = polyfit(x, y, 1);y_fit_linear = polyval(p, x);%多项式曲线拟合p = polyfit(x, y, n);% n 为多项式的次数y_fit_poly = polyval(p, x);%指数曲线拟合p = fit(x, y, 'exp1');y_fit_exp = p(x);```(三)数值积分1、梯形公式```matlabfunction T = trapezoidal_rule(f, a, b, n)h =(b a) / n;x = a:h:b;y = f(x);T = h ((y(1) + y(end))/ 2 + sum(y(2:end 1)));end```2、辛普森公式```matlabfunction S = simpson_rule(f, a, b, n)if mod(n, 2) ~= 0error('n 必须为偶数');endh =(b a) / n;x = a:h:b;y = f(x);S = h / 3 (y(1) + 4 sum(y(2:2:end 1))+ 2 sum(y(3:2:end 2))+ y(end));end```3、柯特斯公式```matlabfunction C = cotes_rule(f, a, b, n)h =(b a) / n;x = a:h:b;y = f(x);w = 7, 32, 12, 32, 7 / 90;C = h sum(w y);end```4、高斯勒让德求积公式```matlabfunction G = gauss_legendre_integration(f, a, b)x, w = gauss_legendre(5);%选择适当的节点数t =(b a) / 2 x +(a + b) / 2;G =(b a) / 2 sum(w f(t));end```(四)数值微分```matlabfunction dydx = numerical_derivative(f, x, h)dydx =(f(x + h) f(x h))/(2 h);end```四、实验结果与分析(一)误差分析通过不同精度的计算,发现随着精度的提高,误差逐渐减小,但计算时间也相应增加。
数值分析在生活中的应用举例及Matlab实现
![数值分析在生活中的应用举例及Matlab实现](https://img.taocdn.com/s3/m/9e9c42dbdb38376baf1ffc4ffe4733687f21fc45.png)
一、最小二乘法,用MATLAB实现1. 数值实例下面给定的是乌鲁木齐最近1个月早晨7:00左右(新疆时间)的天气预报所得到的温度,按照数据找出任意次曲线拟合方程和它的图像。
下面用MATLAB编程对上述数据进行最小二乘拟合。
下面用MATLAB编程对上述数据进行最小二乘拟合2、程序代码x=[1:1:30];y=[9,10,11,12,13,14,13,12,11,9,10,11,12,13,14,12,11,10,9,8,7,8,9,11,9,7,6,5,3,1];a1=polyfit(x,y,3) %三次多项式拟合%a2= polyfit(x,y,9) %九次多项式拟合%a3= polyfit(x,y,15) %十五次多项式拟合%b1=polyval(a1,x)b2=polyval(a2,x)b3=polyval(a3,x)r1= sum((y-b1).^2) %三次多项式误差平方和%r2= sum((y-b2).^2) %九次次多项式误差平方和%r3= sum((y-b3).^2) %十五次多项式误差平方和%plot(x,y,'*') %用*画出x,y图像%hold onplot(x,b1, 'r') %用红色线画出x,b1图像%hold onplot(x,b2, 'g') %用绿色线画出x,b2图像%hold onplot(x,b3, 'b:o') %用蓝色o线画出x,b3图像%3、数值结果不同次数多项式拟合误差平方和为:r1=67.6659r2=20.1060r3=3.7952r1、r2、r3分别表示三次、九次、十五次多项式误差平方和。
4、拟合曲线如下图二、 线性方程组的求解( 高斯-塞德尔迭代算法 )1、实例: 求解线性方程组(见书P233页)⎪⎪⎩⎪⎪⎨⎧=++=-+=+-3612363311420238321321321x x x x x x x x x 记A x=b, 其中⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡=⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡=⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡--=363320,,12361114238321b x A x x x任取初始值()()Tx0000=,进行迭代。
数值研究分析试验幂法与反幂法matlab
![数值研究分析试验幂法与反幂法matlab](https://img.taocdn.com/s3/m/ace113301a37f111f0855b13.png)
数值分析试验幂法与反幂法matlab————————————————————————————————作者:————————————————————————————————日期:一、问题的描述及算法设计(一)问题的描述我所要做的课题是:对称矩阵的条件数的求解设计1、求矩阵A 的二条件数问题 A=⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡----210121012 2、设计内容:1)采用幂法求出A 的. 2)采用反幂法求出A 的.3)计算A 的条件数 ⅡA Ⅱ2* ⅡA -1Ⅱ2=cond2(A )=/.(精度要求为10-6) 3、设计要求1)求出ⅡA Ⅱ2。
2)并进行一定的理论分析。
(二)算法设计1、幂法算法(1)取初始向量u )0((例如取u )0(=(1,1,…1)T ),置精度要求ε,置k=1.(2)计算v )(k =Au )1(-k ,m k =max(v )(k ), u )(k = v )(k / m k(3)若| m k = m 1-k |<ε,则停止计算(m k 作为绝对值最大特征值1λ,u )(k 作为相应的特征向量)否则置k=k+1,转(2)2、反幂法算法(1)取初始向量u )0((例如取u )0(=(1,1,…1)T ),置精度要求ε,置k=1.(2)对A 作LU 分解,即A=LU(3)解线性方程组 Ly )(k =u )1(-k ,Uv )(k =y )(k(4)计算m k =max(v )(k ), u )(k = v )(k / m k(5)若|m k =m 1-k |<ε,则停止计算(1/m k 作为绝对值最小特征值n λ,u )(k 作为相应的特征向量);否则置k=k+1,转(3).二、算法的流程图(一)幂法算法的流程图noyes开始k=0;m1=0 v=A*u [vmax,i]=max(abs(v m=v(i);u=v/mabs(m-m1)< 1e-6 index=1;break; 输出:m ,u ,index 结束m1=m;k=k+1(二)反幂法算法的流程图no yes开始输入A ;[m ,u,index] =pow_inv(A,1e-6) k=0;m1=0 v=invA*u [vmax,i]=max(abs(v)) m=v(i);u=v/mabs(m-m1)< 1e-6 index=1;break; 输出:m ,u ,index结束 m1=m;k=k+1 输入A ;三、算法的理论依据及其推导(一)幂法算法的理论依据及推导幂法是用来确定矩阵的主特征值的一种迭代方法,也即,绝对值最大的特征值。
复合形法matlab程序及例题
![复合形法matlab程序及例题](https://img.taocdn.com/s3/m/0d93af2ec381e53a580216fc700abb68a982ada7.png)
复合形法matlab程序及例题复合形法是一种数值分析方法,用于求解多元函数的最小值。
本文将介绍复合形法的基本原理及其在MATLAB中的实现,同时给出一个具体的例题进行演示。
一、复合形法基本原理复合形法(Nelder-Mead算法)是一种无导数优化方法,它的基本思想是在函数域内随机选择一些点,通过对这些点进行逐步变换、调整,最终找到极小值点。
其算法流程如下:1. 选定n+1个点作为初始点,其中n为被优化函数的自变量个数。
这些点可以随机生成,也可以通过一定的规则生成。
2. 对这些点进行排序,确定最小值点、最大值点和次大值点。
3. 通过一系列变换操作,将最大值点向中间靠拢,同时保持其他点的位置不变。
4. 根据变换后的结果,重新排序,判断是否满足终止条件。
若不满足,继续进行变换操作,直到满足终止条件为止。
终止条件可以是迭代次数、函数值的变化量或某些停止准则的满足等。
二、复合形法MATLAB程序实现复合形法的MATLAB程序可以通过fminsearch函数实现,其语法格式如下:[x,fval,exitflag,output] = fminsearch(fun,x0,options) 其中,fun为被优化的函数,x0为初始点,options为优化选项。
程序返回变量x为最小值点,fval为最小值,exitflag为程序停止状态,output为程序运行信息。
下面是一个简单的例子,用来演示如何使用MATLAB实现复合形法:1. 定义被优化的函数:f(x,y) = (x-2)^2+(y+1)^2+2function z = myfun(x)z = (x(1)-2)^2 + (x(2)+1)^2 + 2;2. 定义初始点和优化选项:x0 = [-1,1];opts = optimset('Display','iter');3. 运行fminsearch函数:[x,fval,exitflag,output] = fminsearch(@myfun,x0,opts);4. 输出结果:x = 1.9999 -0.9999fval = 2.0000exitflag = 1output =iterations: 30funcCount: 58algorithm: 'Nelder-Mead simplex direct search'message: 'Optimization terminated: relative function value decreasing below OPTIONS.TolFun.'最终输出结果表明,通过30次迭代,算法找到了函数的最小值点为x=1.9999,y=-0.9999,函数值为2.0000。
数值分析matlab数值试验
![数值分析matlab数值试验](https://img.taocdn.com/s3/m/c30e0402a6c30c2259019e10.png)
实验一:误差传播及算法稳定性实验1.21、试验程序:function charpt1_2% 误差传播及算法稳定性实验clc;clear all;promps={'请选择递推关系式,若选E1=1/e,En=1-nEn-1,请输入1,若选EN=0,En-1=(1-En)/n,请输入2:'};I=1;while Iresult=inputdlg(promps,'charpt1_2',1,{'1'});Nb=str2num(char(result));if ((Nb~=1)|(Nb~=2))I=0;endend%%%%%%%%%%%%%%%%%I=1;while Iresult=inputdlg('请输入递推步数 n>=1:','charpt1_2',1,{'10'});steps=str2num(char(result));if (steps>0)&(steps==fix(steps)) %% 如果steps大于0且为整数I=0;endend%%%%%%%%%%%%%%%%%result=inputdlg('请输入计算中所采用的有效数字位数n:','charpt1_2',1,{'5'});Sd=str2num(char(result));format long %% 设置显示精度result=zeros(1,steps); %% 存储计算结果err=result; %% 存储计算的绝对误差值func=result; %% 存储用quadl计算的近似值%%%%%%%%%%%%%%%%%%% 用quadl计算积分近似值for n=1:stepsfun=@(x) x.^n.*exp(x-1);func(n)= quadl(fun,0,1);end%%%%%%%%%%%%%%%%%%% 用自定义算法计算if(Nb==1)digits(Sd);result(1)=subs(vpa(1/exp(1)));for n=2:stepsresult(n)=subs(vpa(1-n*result(n-1)));enderr=abs(result-func);elseif(Nb==2)digits(Sd);result(steps)=0;for n=(steps-1):-1:1result(n)=subs(vpa((1-result(n+1))/(n+1)));enderr=abs(result-func);end%%%%%%%%%%%%%%%%%%% 输出结果数值及图像clf;disp('库函数计算值:');disp(sprintf('%e ',func));disp('递推值:');disp(sprintf('%e ',result));disp('误差值:');disp(sprintf('%e ',err));if(Nb==1)plot([1:steps],result,'-rs',[1:steps],func,':k*',[1:steps],err,'-.bo' );elseif(Nb==2)plot([steps:-1:1],result,'-rs',[steps:-1:1],func,':k*',[steps:-1:1],e rr,'-.bo');endxlabel('第n步');ylabel('计算值');legend('自定义算法结果','库函数计算结果','误差值');grid on2、试验结果:选择递推关系式1,递推步数为10,有效数字为5位,计算结果如下:库函数计算值:3.678794e-001 2.642411e-001 2.072766e-001 1.708934e-001 1.455329e-0011.268024e-001 1.123836e-001 1.009323e-001 9.161229e-002 8.387707e-002递推值:3.678800e-001 2.642400e-001 2.072800e-001 1.708800e-001 1.456000e-0011.264000e-001 1.152000e-001 7.840000e-0022.944000e-001 -1.944000e+000误差值:5.588280e-007 1.117662e-006 3.352927e-006 1.341222e-0056.705713e-005 4.023702e-004 2.816427e-003 2.253226e-002 2.027877e-001 2.027877e+00012345678910第n 步计算值选择递推关系式2,递推步数为10,有效数字为5位,计算结果如下: 库函数计算值:3.678794e-001 2.642411e-001 2.072766e-001 1.708934e-001 1.455329e-001 1.268024e-001 1.123836e-001 1.009323e-001 9.161229e-002 8.387707e-002 递推值:3.678800e-001 2.642400e-001 2.072800e-001 1.708900e-001 1.455300e-001 1.267900e-001 1.125000e-001 1.000000e-001 1.000000e-001 0.000000e+000 误差值:5.588280e-007 1.117662e-006 3.352927e-006 3.412224e-006 2.942873e-006 1.237016e-005 1.164270e-004 9.322618e-004 8.387707e-003 8.387707e-002第n 步计算值选择递推关系式1,递推步数为10,有效数字为6位,计算结果如下: 库函数计算值:3.678794e-001 2.642411e-001 2.072766e-001 1.708934e-001 1.455329e-001 1.268024e-001 1.123836e-001 1.009323e-001 9.161229e-002 8.387707e-002 递推值:3.678790e-001 2.642420e-001 2.072740e-001 1.709040e-001 1.454800e-001 1.271200e-001 1.101600e-001 1.187200e-001 -6.848000e-002 1.684800e+000 误差值:4.411720e-007 8.823378e-007 2.647073e-006 1.058778e-0055.294287e-005 3.176298e-004 2.223573e-003 1.778774e-002 1.600923e-001 1.600923e+00012345678910第n 步计算值选择递推关系式2,递推步数为10,有效数字为6位,计算结果如下: 库函数计算值:3.678794e-001 2.642411e-001 2.072766e-001 1.708934e-001 1.455329e-001 1.268024e-001 1.123836e-001 1.009323e-001 9.161229e-002 8.387707e-002 递推值:3.678800e-001 2.642410e-001 2.072770e-001 1.708930e-001 1.455360e-001 1.267860e-001 1.125000e-001 1.000000e-001 1.000000e-001 0.000000e+000 误差值:5.588280e-007 1.176622e-007 3.529274e-007 4.122239e-007 3.057127e-006 1.637016e-005 1.164270e-004 9.322618e-004 8.387707e-003 8.387707e-002第n 步计算值3、结果分析:很明显第二种递推式结果要比第一种好,式1在第七步后有明显误差,而式2在第三步后基本与近似解一致。
matlab幂法求特征值和特征向量方法实现和函数表示
![matlab幂法求特征值和特征向量方法实现和函数表示](https://img.taocdn.com/s3/m/5bd971af112de2bd960590c69ec3d5bbfc0ada7e.png)
matlab幂法求特征值和特征向量方法实现和函数表示1. 引言在数值分析中,求解特征值和特征向量是一项重要而且经常出现的任务。
特征值和特征向量在矩阵和线性代数中有着广泛的应用,涉及到许多领域,如机器学习、信号处理、结构动力学等。
在matlab中,幂法是一种常用的求解特征值和特征向量的方法,同时也有对应的函数可以实现这一过程。
2. 幂法的原理幂法是一种迭代方法,它利用矩阵的特征值和特征向量的性质,通过不断地迭代计算,逼近矩阵的主特征值和对应的特征向量。
具体来说,假设A是一个n阶矩阵,它的特征值λ1>λ2≥...≥λn,并且对应着线性无关的特征向量v1,v2,...,vn。
如果选择一个任意的非零初始向量x0,并进行以下迭代计算:```x(k+1) = Ax(k) / ||Ax(k)||```其中,||.||表示向量的模长。
不断迭代计算后,x(k)将收敛到矩阵A的主特征向量v1上,并且相应的特征值即为A的主特征值λ1。
3. matlab实现幂法求解特征值和特征向量在matlab中,幂法的实现也非常简单。
可以使用自带的eig函数,该函数可以直接求解矩阵的特征值和特征向量。
使用方法如下:```[V,D] = eig(A)```其中,A为待求解的矩阵,V为特征向量矩阵,D为特征值矩阵。
利用eig函数,即可一步到位地求解矩阵的特征值和特征向量,非常简单方便。
4. 函数表示幂法求解特征值和特征向量的过程可以表示为一个matlab函数。
通过封装相关的迭代算法和收敛判据,可以方便地实现幂法的函数表示。
可以定义一个名为powerMethod的函数:```matlabfunction [lambda, v] = powerMethod(A, x0, maxIter, tol)% 初始化k = 1;x = x0;% 迭代计算while k <= maxItery = A * x;lambda = norm(y, inf);x = y / lambda;% 检查收敛性if norm(A * x - lambda * x) < tolbreak;endk = k + 1;endv = x;end```利用这个函数,就可以自己实现幂法求解特征值和特征向量的过程。
matlab数值分析实验报告
![matlab数值分析实验报告](https://img.taocdn.com/s3/m/f7540277effdc8d376eeaeaad1f34693daef10e2.png)
matlab数值分析实验报告Matlab数值分析实验报告引言数值分析是一门研究利用计算机进行数值计算和模拟的学科,它在科学计算、工程技术和金融等领域有着广泛的应用。
本次实验报告将介绍在Matlab环境下进行的数值分析实验,包括数值微分、数值积分和线性方程组求解等内容。
一、数值微分数值微分是通过数值方法计算函数的导数,常用的数值微分方法有前向差分、后向差分和中心差分。
在Matlab中,可以使用diff函数来计算函数的导数。
例如,对于函数f(x)=x^2,在Matlab中可以使用如下代码进行数值微分的计算:```matlabsyms x;f = x^2;df = diff(f, x);```二、数值积分数值积分是通过数值方法计算函数的定积分,常用的数值积分方法有梯形法则、辛普森法则和龙贝格积分法。
在Matlab中,可以使用trapz、quad和integral等函数来进行数值积分的计算。
例如,对于函数f(x)=sin(x),可以使用如下代码进行数值积分的计算:```matlabx = linspace(0, pi, 100);y = sin(x);integral_value = trapz(x, y);```三、线性方程组求解线性方程组求解是数值分析中的重要问题,常用的求解方法有高斯消元法和LU 分解法。
在Matlab中,可以使用\操作符来求解线性方程组。
例如,对于线性方程组Ax=b,可以使用如下代码进行求解:```matlabA = [1, 2; 3, 4];b = [5; 6];x = A\b;```四、实验结果与分析在本次实验中,我们分别使用Matlab进行了数值微分、数值积分和线性方程组求解的计算。
通过实验结果可以发现,Matlab提供了丰富的数值计算函数和工具,能够方便地进行数值分析的计算和求解。
数值微分的计算结果与解析解相比较,可以发现数值微分的误差随着步长的减小而减小,但是当步长过小时,数值微分的误差会受到舍入误差的影响。
MATLAB 数值分析
![MATLAB 数值分析](https://img.taocdn.com/s3/m/ea622cc1aa00b52acfc7cab1.png)
13.2
极小化
作图除了提供视觉信息外,还常常需要确定一个函数的其它更多的特殊属性。在许多 应用中,特别感兴趣的是确定函数的极值,即最大值(峰值)和最小值(谷值)。数学上, 可通过确定函数导数(斜率)为零的点,解析上求出这些极值点。检验 humps 的图形在峰 值和谷值点上的斜率就很容易理解这个事实。显然,如果定义的函数简单,则这种方法常 常奏效。然而,即使很多容易求导的函数,也常常很难找到导数为零的点。在这种情况下, 以及很难或不可能解析上求得导数的情况下,必须数值上寻找函数的极值点。MATLAB 提 供了两个完成此功能的函数 fmin 和 fmins。 这两个函数分别寻找一维或 n 维函数的最小值。 这里仅讨论 fmin。有关 fmins 的详细信息,参阅《MATLAB 参考指南》。因为 f(x)的最大 值等于-f(x)的最小值,所以,上述 fmin 和 fmins 可用来求最大值和最小值。如果还不清楚, 把上述图形倒过来看,在这个状态下,峰值变成了谷值,而谷值则变成了峰值。 为了解释求解一维函数的最小值和最大值, 再考虑上述例子。 从图 13.2 可知, 在 xmax=0.7 附 近 有 一 个 最 大 值 , 并 且 在 xmin=4 附 近 有 一 个 最 小 值 。 而 这 些 点 的 解 析 值 为 : x m a x / 4 0.785 和 x min 5 / 4 393 . 。为了方便,用文本编辑器编写一个脚本 M 文件,并用 fmin 寻出数值上极值点,给出函数主体如下: % ex_fmin.m fn=‘ 2*exp(-x)*sin(x) ‘; xmin=fmin(fn , 2 , 5)
% better approximation
自然地,上述两个结果不同。基于对图形的观察,粗略近似可能低估了实际面积。除 非特别精确,没有准则说明哪种近似效果更好。很明显,如果人们能够以某种方式改变单 个梯形的宽度,以适应函数的特性,即当函数变化快时,使得梯形的宽度变窄,这样就能 够得到更精确的结果。 MATLAB 的函数 quad 和 quad8 是基于数学上的正方形概念来计算函数的面积, 这些 积分函数的操作方式一样。为获得更准确的结果,两个函数在所需的区间都要计算被积函 数。 此外, 与简单的梯形比较, 这两个函数进行更高阶的近似, 而且 quad8 比 quad 更精确。 这两个函数的调用方法与 fzero 相同,即 >>area=quad(‘ humps ‘ , -1 , 2) % find area between -1 and 2
matlab 数值解
![matlab 数值解](https://img.taocdn.com/s3/m/94c62b0111661ed9ad51f01dc281e53a5902514e.png)
matlab 数值解Matlab 数值解Matlab 是一种强大的数学软件,它包含了很多数学工具箱,可以用于数值分析和求解数学问题。
在本文中,我们将介绍Matlab 中的数值解方法,包括数值积分、数值微分、非线性方程求解和常微分方程的数值解法。
数值积分数值积分是一种数学方法,用于求解函数的定积分。
在Matlab 中,可以使用 quad 和 quadl 函数进行数值积分。
其中,quad 函数用于计算一般积分,而 quadl 函数用于计算不定积分。
数值微分数值微分是一种数学方法,用于计算函数的导数。
在Matlab 中,可以使用diff 和gradient 函数进行数值微分。
其中,diff 函数用于计算一维函数的导数,而 gradient 函数用于计算多维函数的梯度。
非线性方程求解非线性方程是一种形式为 f(x)=0 的方程,其中 f(x) 是一个非线性函数。
在 Matlab 中,可以使用 fzero 和 fsolve 函数进行非线性方程求解。
其中,fzero 函数用于求解单变量非线性方程,而fsolve 函数用于求解多变量非线性方程。
常微分方程的数值解法常微分方程是一种形式为y'=f(t,y) 的方程,其中y 是未知函数,t 是自变量,f(t,y) 是已知函数。
在Matlab 中,可以使用ode45 和ode23 函数进行常微分方程的数值解法。
其中,ode45 函数是一种常用的数值解法,可以求解大部分常微分方程,而 ode23 函数则是一种高效的数值解法,适用于求解简单的常微分方程。
总结在本文中,我们介绍了Matlab 中的数值解方法,包括数值积分、数值微分、非线性方程求解和常微分方程的数值解法。
这些方法可以帮助我们快速、准确地求解数学问题,提高数学建模的效率和精度。
同济大学数值分析matlab编程
![同济大学数值分析matlab编程](https://img.taocdn.com/s3/m/f6f8c9fc5ff7ba0d4a7302768e9951e79b89692f.png)
同济⼤学数值分析matlab编程MATLAB 编程题库1.下⾯的数据表近似地满⾜函数21cx bax y ++=,请适当变换成为线性最⼩⼆乘问题,编程求最好的系数c b a ,,,并在同⼀个图上画出所有数据和函数图像.625.0718.0801.0823.0802.0687.0606.0356.0995.0628.0544.0008.0213.0362.0586.0931.0ii y x ----解:>> x=[-0.931 -0.586 -0.362 -0.213 0.008 0.544 0.628 0.995]'; >> y=[0.356 0.606 0.687 0.802 0.823 0.801 0.718 0.625]'; >> A= [x ones(8,1) -x.^2.*y]; >> z=A\y;>> a=z(1); b=z(2); c=z(3); >>xh=[-1:0.1:1]';>>yh=(a.*xh+b)./(1+c.*xh.^2); >>plot(x,y,'r+',xh,yh,'b*')2.若在Matlab ⼯作⽬录下已经有如下两个函数⽂件,写⼀个割线法程序,求出这两个函数精度为1010-的近似根,并写出调⽤⽅式:>> edit gexianfa.mfunction [x iter]=gexianfa(f,x0,x1,tol) iter=0;x=x1;while(abs(feval(f,x))>tol) iter=iter+1;x=x1-feval(f,x1).*(x1-x0)./(feval(f,x1)-feval(f,x0)); x0=x1;x1=x; end>> edit f.m function v=f(x) v=x.*log(x)-1;>> edit g.m function z=g(y) z=y.^5+y-1;>> [x1 iter1]=gexianfa('f',1,3,1e-10) x1 =1.7632 iter1 = 6>> [x2 iter2]=gexianfa('g',0,1,1e-10) x2 =0.7549 iter2 = 83.使⽤GS 迭代求解下述线性代数⽅程组:123123123521242103103x x x x x x x x x ì++=--++=í???-+=??解:>> edit gsdiedai.mfunction [x iter]=gsdiedai(A,x0,b,tol) D=diag(diag(A)); L=D-tril(A); U=D-triu(A); iter=0; x=x0;>> A=[5 2 1;-1 4 2;1 -3 10]; >> b=[-12 10 3]'; >>tol=1e-4; >>x0=[0 0 0]';>> [x iter]=gsdiedai(A,x0,b,tol); >>x x =-3.0910 1.2372 0.9802 >>iter iter = 64.⽤四阶Range-kutta ⽅法求解下述常微分⽅程初值问题(取步长h=0.01),(1)2x dy y e xy dx y ì??=++?í??=??解:>> edit ksf2.mfunction v=ksf2(x,y)v=y+exp(x)+x.*y; >> a=1;b=2;h=0.01; >> n=(b-a)./h; >> x=[1:0.01:2]; >>y(1)=2;>>for i=2:(n+1)k1=h*ksf2(x(i-1),y(i-1));k2=h*ksf2(x(i-1)+0.5*h,y(i-1)+0.5*k1); k3=h*ksf2(x(i-1)+0.5*h,y(i-1)+0.5*k2); k4=h*ksf2(x(i-1)+h,y(i-1)+k3); y(i)=y(i-1)+(k1+2*k2+2*k3+k4)./6; end >>y调⽤函数⽅法>> edit Rangekutta.mfunction [x y]=Rangekutta(f,a,b,h,y0) x=[a:h:b]; n=(b-a)/h; y(1)=y0; for i=2:(n+1)k1=h*(feval(f,x(i-1),y(i-1)));k2=h*(feval(f,x(i-1)+0.5*h,y(i-1)+0.5*k1)); k3=h*(feval(f,x(i-1)+0.5*h,y(i-1)+0.5*k2)); k4=h*(feval(f,x(i-1)+h,y(i-1)+k3)); y(i)=y(i-1)+ (k1+2*k2+2*k3+k4)./6; end>> [x y]=Rangekutta('ksf2',1,2,0.01,2); >>y5.取0.2h =,请编写Matlab 程序,分别⽤欧拉⽅法、改进欧拉⽅法在12x ≤≤上求解初值问题。
数值分析matlab方法插值法
![数值分析matlab方法插值法](https://img.taocdn.com/s3/m/20873384af45b307e9719791.png)
其中,
n
【注】
x [a, b] w(x) (x x j ) j 0
(1)误差估计
Rn (x)
M n1 (n 1)!
w( x)
M n1
max
x( a ,b )
f
(n1) (x)
(2)余项与 x、M n1 节点的位置、个数 n 有关
(3)当 f (x)是 n 的多项式时Ln (x) f (x) n
M2 2!
(x
x0 )(x
x1 )
其中,
M2
max
x( x0 , x1 )
f
(x)
x
[
6
,
4
]
,所以
R1
(
5
24
)
sin
2!
4 (5
24
)( 5
6 24
)
4
0.0061
2)
抛物插值误差估计.因为
R2 (x)
M3 3!
(x
x0 )(x
x1)(x
x2 )
其中,
M3
max
x( x0 ,x2 )
f (x)
yiynewtonbackwardxyxicos035yi09394数值分析插值法55埃尔米特插值2n12n2数值分析插值法551埃尔米特插值多项式的存在唯一性2n1数值分析插值法数值分析插值法552埃尔米特插值余项553三次埃尔米特插值多项式maxsinxsin1数值分析插值法56561高次插值的病态性质0908原函数150706050405分段线性插值0302543214321055数值分析插值法562分段低次插值方法563分段低次插值余项090807060504030201090807060504030201543214321数值分析插值法57三次样条插值571三次样条插值572三弯矩法数值分析插值法573三次样条插值的误差估计与收敛性58插值运算的matlab函数581一维插值函数interp1yiinterp1xyximethod?linear?yiinterp1xyxilinear?1200135019
数值分析(Hilbert矩阵)病态线性方程组的求解Matlab程序
![数值分析(Hilbert矩阵)病态线性方程组的求解Matlab程序](https://img.taocdn.com/s3/m/4d383598690203d8ce2f0066f5335a8102d266a1.png)
数值分析(Hilbert矩阵)病态线性⽅程组的求解Matlab程序(Hilbert 矩阵)病态线性⽅程组的求解理论分析表明,数值求解病态线性⽅程组很困难。
考虑求解如下的线性⽅程组的求解Hx = b ,期中H 是Hilbert 矩阵,()ij n n H h ?=,11ij h i j = +-,i ,j = 1,2,…,n 1. 估计矩阵的2条件数和阶数的关系2. 对不同的n ,取(1,1,,1)n x =∈,分别⽤Gauss 消去,Jacobi 迭代,Gauss-seidel 迭代,SOR 迭代和共轭梯度法求解,⽐较结果。
3. 结合计算结果,试讨论病态线性⽅程组的求解。
第1⼩题:condition.m %第1⼩题程序t1=20;%阶数n=20x1=1:t1;y1=1:t1;for i=1:t1H=hilb(i);y1(i)=log(cond(H));endplot(x1,y1);xlabel('阶数n');ylabel('2-条件数的对数(log(cond(H))');title('2-条件数的对数(log(cond(H))与阶数n 的关系图');t2=200;%阶数n=200x2=1:t2;y2=1:t2;for i=1:t2H=hilb(i);y2(i)=log(cond(H));endplot(x2,y2);xlabel('阶数n');ylabel('2-条件数的对数(log(cond(H))');title('2-条件数的对数(log(cond(H))与阶数n 的关系图');画出Hilbert 矩阵2-条件数的对数和阶数的关系n=200时n=20时从图中可以看出,1)在n⼩于等于13之前,图像近似直线log(cond(H))~1.519n-1.8332)在n⼤于13之后,图像趋于平缓,并在⼀定范围内上下波动,同时随着n的增加稍有上升的趋势第2⼩题:solve.m%m第2⼩题主程序N=4000;。
- 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。
数值分析matlab实现高斯消元法:function[RA,RB,n,X]=gaus(A,b)B=[A b];n=length(b);RA=rank(A);RB=rank(B);zhica=RB-RA;if zhica>0,disp('请注意:因为RA~=RB,所以此方程组无解.')returnendif RA==RBif RA==ndisp('请注意:因为RA=RB=n,所以此方程组有唯一解.')X=zeros(n,1);C=zeros(1,n+1);for p=1:n-1for k=p+1:nm=B(k,p)/B(p,p);B(k,p:n+1)=B(k,p:n+1)-m*B(p,p:n+1); endendb=B(1:n,n+1);A=B(1:n,1:n);X(n)=b(n)/A(n,n);for q=n-1:-1:1X(q)=(b(q)-sum(A(q,q+1:n)*X(q+1:n)))/A(q,q);endelsedisp('请注意:因为RA=RB<n,所以此方程组有无穷多解.')endendX列主元消去法:function[RA,RB,n,X]=liezhu(A,b)B=[A b];n=length(b);RA=rank(A);RB=rank(B);zhica=RB-RA;if zhica>0,disp('请注意:因为RA~=RB,所以此方程组无解.')returnendif RA==RBif RA==ndisp('请注意:因为RA=RB=n,所以此方程组有唯一解.')X=zeros(n,1);C=zeros(1,n+1);for p=1:n-1[Y,j]=max(abs(B(p:n,p)));C=B(p,:);B(p,:)=B(j+p-1,:);B(j+p-1,:)=C;for k=p+1:nm=B(k,p)/B(p,p);B(k,p:n+1)=B(k,p:n+1)-m*B(p,p:n+1);endendb=B(1:n,n+1);A=B(1:n,1:n);X(n)=b(n)/A(n,n);for q=n-1:-1:1X(q)=(b(q)-sum(A(q,q+1:n)*X(q+1:n)))/A(q,q);endelsedisp('请注意:因为RA=RB<n,所以此方程组有无穷多解.')endendXJacobi迭代法:例1用jacobi迭代法求解代数线性代数方程组,保留四位有效数字(err=1e-4)其中A=[8-11;2 10-1;11-5];b=[1;4;3]。
解:编写jacobi迭代法的函数文件,保存为jacobi.mfunction[x,k]=jacobi(A,b,x0,eps,N)%求解Ax=b;x0为初始列向量;eps为误差容限;N为最大迭代次数%输出x为近似解;k为迭代次数n=length(A);x=zeros(n,1);for k=1:Nfor i=1:n―――――――endif norm(x-x0,inf)<eps%迭代终止条件%if(max(abs(x-x0)))<epsbreak;endx0=x;end编写主程序如下format longclearA=[8-11;210-1;11-5];b=[1;4;3];x0=[0.125;0.4;-0.6];%x0为初始列向量N为最大迭代次数err=1e-4;%err为误差容限N=25;%N为最大迭代次数[x,k]=jacobi(A,b,x0,err,N)得到结果如下x=0.224923156250000.30561995000000-0.49388680000000k=6024681012HourD e g r e e s C e l s i u sGauss-seidel 迭代法:例2用Gauss-seidel 迭代法求解代数线性代数方程组,保留四位有效数字(err =1e-4)其中A=[8-11;210-1;11-5];b=[1;4;3]。
解:编写Gauss-seidel 迭代法的函数文件,保存为gaus.mfunction [x,k]=gaus(A,b,x0,eps,N)%求解Ax=b ;x0为初始列向量;eps 为误差容限;N 为最大迭代次数%输出x 为近似解;k 为迭代次数n=length(A);x=zeros(n,1);for k=1:Nfor i=1:n――――――endif norm(x-x0,inf)<eps %迭代终止条件%if (max(abs(x-x0)))<epsbreak;end x0=x;end编写主程序如下format long clearA=[8-11;210-1;11-5];b=[1;4;3];x0=[0.125;0.4;-0.6];%x0为初始列向量N 为最大迭代次数err=1e-4;%err 为误差容限N=25;%N 为最大迭代次数[x,k]=gaus(A,b,x0,err,N)输出结果为x =0.224939378906250.30562326171875-0.49388747187500k =5用matlab 做插值计算一维插值函数:Yi=interp1(x,y,xi,’method’)其中yi ——xi 处的插值结果;x,y ——插值节点;xi ——被插值点;’method’——插值方法,分为’nearest’—最邻近插值;’linear’—线性插值;’spline’—三次样条插值;’cubic’—立方插值。
例如:>>hours=1:12;>>temps=[589152529313022252724];>>h=1:0.1:12;>>t=interp1(hours,temps,h,'linear');>>plot(hours,temps)>>xlabel('Hour'),ylabel('Degrees Celsius')二维插值:1,网格节点数据的插值24681012z=interp2(x0,y0,z0,x,y,’method’)其中,z ——被插值点的函数值;x0,y0,z0——插值节点;x,y ——被插值点;’method’——插值方法,与一维插值相同。
2,散点数据的插值计算cz=griddata(x,y,z,cx,cy,’method’)与网格节点数据相似。
Lagrange 多项式插值(n 次):function f=Lagrange(x,fx,inx)n=length(x);m=length(inx);for i=1:m;z=inx(i);s=0.0;for k=1:np=1.0;for j=1:nif j~=kp=p*(z-x(j))/(x(k)-x(j));end end s=p*fx(k)+s;end f(i)=s;endplot(x,fx,'O',inx,f)例:例:已知数据如下,求1—12之间以0.2为间隔的所有值,并画出图形:xi=123456789101112;fxi=1223434-13425233494523。
解:图像如右图。
>>x=[1:12];>>fx=[1223434-13425233494523];>>xi=[1:0.2:12];>>Lagrange(x,fx,xi)ans =Columns 1through 912.0000-60.593718.2765124.9778202.5952234.0000223.3757184.1249131.4738Columns 10through 1878.425334.0000 2.9467-13.6885-17.5810-12.0379-1.000011.755623.1624Columns 19through 2731.161134.773034.000029.605422.833215.11537.8099 2.0000-1.6307Columns 28through 36-2.8397-1.7907 1.0404 5.00009.402413.664317.403320.483423.0000Columns37through4525.203727.376929.685832.040034.000034.774233.342628.7320 20.4439Columns46through549.0000-3.4848-12.8605-12.8873 4.059245.0000112.3788197.1817 267.9699Columns55through56254.343923.0000Hermite插值法:function f=Hermite(x,y,y_1,x0)syms t;f=0.0;if(length(x)==length(y))if(length(y)==length(y_1))n=length(x);elsedisp('y和y的导数的维数不相等!');return;endelsedisp('x和y的维数不相等!');return;endfor i=1:nh=1.0;a=0.0;for j=1:nif(j~=i)h=h*(t-x(j))^2/((x(i)-x(j))^2);a=a+1/(x(i)-x(j));endendf=f+h*((x(i)-t)*(2*a*y(i)-y_1(i))+y(i));if(i==n)if(nargin==4)f=subs(f,'t',x0);elsef=vpa(f,6);endendend三次样条插值:调用Matlab自带函数:spline。