导航:首页 > 源码编译 > matlab拉格朗日插值算法

matlab拉格朗日插值算法

发布时间:2022-09-13 16:23:24

① 拉格朗日插值的matlab代码

1、给出一列数据之后,作图如下:aa= randn(100,1);plot(aa);。

② 在MATLAB编程实验中,用拉格朗日插值法跟牛顿插值法运行之后计算的结果为什么是一样的

根据插值多项式的唯一性,两种方法的结果应该是一样的。条条道路通罗马,只是方法不同而已,牛顿法要比拉格朗日法优越简单。

Matlab函数M文件Lagrange程序function yy=lagrange(x,y,xi) m=length(x)上面是拉格朗日插值法,其中xi为要计算的数值比如 x=[0 3 5 9 31];Q

clear all;clc

x0=1:5;

y0=sin(x0);

x=1:0.2:2;

y0=lagrange(x0,y0,x)

命令窗口输这个就没有问题。

(2)matlab拉格朗日插值算法扩展阅读:

如果这特定函数是多项式,就称它为插值多项式。利用插值基函数很容易得到拉格朗日插值多项式,公式结构紧凑,在理论分析中甚为方便,但当插值节点增减时全部插值基函数均要随之变化,整个公式也将发生变化,这在实际计算中是很不方便的,为了克服这一缺点,提出了牛顿插值。

③ 怎么用matlab利用拉格朗日插值计算法的原理编写并计算函数所在节点的近似值。 谢谢。

.m文件
function yy=lagrange(x1,y1,xx)
%本程序为Lagrange1插值,其中x1,y1
%为插值节点和节点上的函数值,输出为插值点xx的函数值,
%xx可以是向量。
syms x
n=length(x1);
for i=1:n
t=x1;t(i)=[];L(i)=prod((x-t)./(x1(i)-t));% L向量用来存放插值基函数
end
u=sum(L.*y1);
p=simplify(u) % p是简化后的Lagrange插值函数(字符串)
yy=subs(p,x,xx);
clf
plot(x1,y1,'ro',xx,yy,'*')
====================================
x=[ 0.4 0.5 0.6 0.7 0.8];
y=[-0.916291;-0.693147;-0.510826;-0.356675;-0.223144]';
yy=lagrange(x,y,0.54)

p =

- (14363668061545223*x^4)/6755399441055744 + (229230406283396627*x^3)/33776997205278720 - (6086876668119665137*x^2)/675539944105574400 + (23595121244981107513*x)/3377699720527872000 - 186390055565518223/70368744177664000

yy =

-0.6161

④ MATLAB中如何利用拉格朗日插值法作图

计算卫星定位么?? 要几次插值的 给你一次和二次的吧
1.n个节点分段Lagrange插值多项式;
%2.使用格式y=lagrange(x0,y0,x,k);
%3.输入项x0为n维插值节点向量,y0为n维被插函数值向量;
%4.x为m维插值点向量,k为分段插值多项式次数,不超过3,缺省为k=1;
%5.输出y为插值点x处的函数值;
%6.本程序于2002.4.21.编写?
function y=lagrange(x0,y0,x,k)
if nargin<4,k=1;end
if k>3,error('分段次数过高,容易产生Runge现象,请重新选择次数k.'),end
n=length(x0);
m=length(x);
nn=1;
for i=1:m
u=x(i);
switch k
%---------------------------针对不同的k判断插值区间
case 1
if u<=x0(2)
t=1;
elseif u>x0(n-1)
t=n-1;
else
for j=nn+1:n-2
if u>x0(j)&u<=x0(j+1)
t=j;nn=t-1;break
end
end
end
%---------------------------
case 2
if u<=(x0(2)+x0(3))/2
t=1;
elseif u>(x0(n-2)+x0(n-1))/2
t=n-2;
else
for j=nn+1:n-3
if u>(x0(j)+x0(j+1))/2 & u<=(x0(j+1)+x0(j+2))/2
t=j;nn=t-1;break
end
end
end
%---------------------------
case 3
if u<=x0(3)
t=1;
elseif u>x0(n-2)
t=n-3;
else
for j=nn+2:n-3
if u>x0(j) & u<=x0(j+1)
t=j-1;nn=t;break
end
end
end
end
%--------------------------------主程序
s=0;
for j=t:t+k
L=1;
for p=t:t+k
if p~=j
L=L*(u-x0(p))/(x0(j)-x0(p));
end
end
s=s+L*y0(j);
end
y(i)=s;
end

分段线性Lagrange插值
% 命令格式:y=lagrange1(x0,y0,x)
% x0为节点向量,y0为对应的函数值向量,
% x为插值点向量,返回值y为x处的函数近似值向量。
function y=lagrange1(x0,y0,x)
[n1,n]=size(x0);[n1,m]=size(x);
for i=1:m
u=x(i);
if u<=x0(2)
y(i)=y0(1)*(u-x0(2))/(x0(1)-x0(2))+...
y0(2)*(u-x0(1))/(x0(2)-x0(1));
elseif u>=x0(n-1)
y(i)=y0(n-1)*(u-x0(n))/(x0(n-1)-x0(n))+...
y0(n)*(u-x0(n-1))/(x0(n)-x0(n-1));
else
for k=2:n-1
if u>=x0(k) & u<=x0(k+1)
y(i)=y0(k)*(u-x0(k+1))/(x0(k)-x0(k+1))+...
y0(k+1)*(u-x0(k))/(x0(k+1)-x0(k));
end
end
end
end
y;

⑤ 大神 求解 拉格朗日插值 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

SOR迭代法的Matlab程序

function [x]=SOR_iterative(A,b)
% 用SOR迭代求解线性方程组,矩阵A是方阵
x0=zeros(1,length(b)); % 赋初值
tol=10^(-2); % 给定误差界
N=1000; % 给定最大迭代次数
[n,n]=size(A); % 确定矩阵A的阶
w=1; % 给定松弛因子
k=1;
% 迭代过程
while k<=N
x(1)=(b(1)-A(1,2:n)*x0(2:n)')/A(1,1);
for i=2:n
x(i)=(1-w)*x0(i)+w*(b(i)-A(i,1:i-1)*x(1:i-1)'-A(i,i+1:n)*x0(i+1:n)')/A(i,i);
end
if max(abs(x-x0))<=tol
fid = fopen('SOR_iter_result.txt', 'wt');
fprintf(fid,'\n********用SOR迭代求解线性方程组的输出结果********\n\n');
fprintf(fid,'迭代次数: %d次\n\n',k);
fprintf(fid,'x的值\n\n');
fprintf(fid, '%12.8f \n', x);
break;
end
k=k+1;
x0=x;
end
if k==N+1
fid = fopen('SOR_iter_result.txt', 'wt');
fprintf(fid,'\n********用SOR迭代求解线性方程组的输出结果********\n\n');
fprintf(fid,'迭代次数: %d次\n\n',k);
fprintf(fid,'超过最大迭代次数,求解失败!');
fclose(fid);
end
Matlab中龙格-库塔(Runge-Kutta)方法原理及实现龙格-库塔(Runge-Kutta)方法是一种在工程上应用广泛的高精度单步算法。由于此算法精度高,采取措施对误差进行抑制,所以其实现原理也较复杂。该算法是构建在数学支持的基础之上的。龙格库塔方法的理论基础来源于泰勒公式和使用斜率近似表达微分,它在积分区间多预计算出几个点的斜率,然后进行加权平均,用做下一点的依据,从而构造出了精度更高的数值积分计算方法。如果预先求两个点的斜率就是二阶龙格库塔法,如果预先取四个点就是四阶龙格库塔法。一阶常微分方程可以写作:y'=f(x,y),使用差分概念。
(Yn+1-Yn)/h= f(Xn,Yn)推出(近似等于,极限为Yn')
Yn+1=Yn+h*f(Xn,Yn)
另外根据微分中值定理,存在0<t<1,使得
Yn+1=Yn+h*f(Xn+th,Y(Xn+th))
这里K=f(Xn+th,Y(Xn+th))称为平均斜率,龙格库塔方法就是求得K的一种算法。
利用这样的原理,经过复杂的数学推导(过于繁琐省略),可以得出截断误差为O(h^5)的四阶龙格库塔公式:
K1=f(Xn,Yn);
K2=f(Xn+h/2,Yn+(h/2)*K1);
K3=f(Xn+h/2,Yn+(h/2)*K2);
K4=f(Xn+h,Yn+h*K3);
Yn+1=Yn+h*(K1+2K2+2K3+K4)*(1/6);
所以,为了更好更准确地把握时间关系,应自己在理解龙格库塔原理的基础上,编写定步长的龙格库塔函数,经过学习其原理,已经完成了一维的龙格库塔函数。
仔细思考之后,发现其实如果是需要解多个微分方程组,可以想象成多个微分方程并行进行求解,时间,步长都是共同的,首先把预定的初始值给每个微分方程的第一步,然后每走一步,对多个微分方程共同求解。想通之后发现,整个过程其实很直观,只是不停的逼近计算罢了。编写的定步长的龙格库塔计算函数:
function [x,y]=runge_kutta1(ufunc,y0,h,a,b)%参数表顺序依次是微分方程组的函数名称,初始值向量,步长,时间起点,时间终点(参数形式参考了ode45函数)
n=floor((b-a)/h);%求步数
x(1)=a;%时间起点
y(:,1)=y0;%赋初值,可以是向量,但是要注意维数
for ii=1:n

x(ii+1)=x(ii)+h;

k1=ufunc(x(ii),y(:,ii));

k2=ufunc(x(ii)+h/2,y(:,ii)+h*k1/2);

k3=ufunc(x(ii)+h/2,y(:,ii)+h*k2/2);

k4=ufunc(x(ii)+h,y(:,ii)+h*k3);

y(:,ii+1)=y(:,ii)+h*(k1+2*k2+2*k3+k4)/6;
%按照龙格库塔方法进行数值求解
end
调用的子函数以及其调用语句:
function dy=test_fun(x,y)
dy = zeros(3,1);%初始化列向量
dy(1) = y(2) * y(3);
dy(2) = -y(1) + y(3);
dy(3) = -0.51 * y(1) * y(2);
对该微分方程组用ode45和自编的龙格库塔函数进行比较,调用如下:
[T,F] = ode45(@test_fun,[0 15],[1 1 3]);
subplot(121)
plot(T,F)%Matlab自带的ode45函数效果
title('ode45函数效果')
[T1,F1]=runge_kutta1(@test_fun,[1 1 3],0.25,0,15);%测试时改变test_fun的函数维数,别忘记改变初始值的维数
subplot(122)
plot(T1,F1)%自编的龙格库塔函数效果
title('自编的 龙格库塔函数')

⑥ 拉格朗日插值 MATLAB

额,你这是lagrange插值?
这里有个lagrange插值的例子,你参照一下吧。
function y=lagrange_interp(xdata,ydata,x)
% Lagrange插值
% 输入参数:
% ---xdata:给定的节点横坐标
% ---ydata:给定的节点纵坐标
% ---x:需要进行插值的节点横坐标
% 输出参数:
% ---y:Lagrange插值函数在x处的函数值
n=length(xdata);m=length(ydata);
if n~=m
error('插值数据长度不等!');
end
ii=1:n;y=zeros(size(x));
for i=ii
ij=find(ii~=i);V=1;
for j=1:length(ij)
if abs(xdata(i)-xdata(ij(j)))<eps
error('输入的n+1个节点不是互异的。');
end
V=V.*(x-xdata(ij(j)));
end
y=y+V*ydata(i)/prod(xdata(i)-xdata(ij));
end

⑦ 用matlab编写拉格朗日插值算法的程序

做了一个测试,希望有所帮助。代码:% 用matlab编写拉格朗日插值算法的程序,并以下面给出的函数表为数据基础,
% 在整个插值区间上采用拉格朗日插值法计算f(0.6),写出程序源代码,输出计算结果
% x -2.15 -1.00 0.01 1.02 2.03 3.25
% y 17.03 7.24 1.05 2.03 17.06 23.05
function main()
clc;
x = [-2.15 -1.00 0.01 1.02 2.03 3.25];
y = [17.03 7.24 1.05 2.03 17.06 23.05 ];
x0 = 0.6;
f = Language(x,y,x0)function f = Language(x,y,x0)
%求已知数据点的拉格朗日插值多项式
%已知数据点的x坐标向量: x
%已知数据点的y坐标向量: y
%插值点的x坐标: x0
%求得的拉格朗日插值多项式或在x0处的插值: fsyms t l;
if(length(x) == length(y))
n = length(x);
else
disp('x和y的维数不相等!');
return; %检错
endh=sym(0);
for (i=1:n)
l=sym(y(i));
for(j=1:i-1)
l=l*(t-x(j))/(x(i)-x(j));
end;
for(j=i+1:n)
l=l*(t-x(j))/(x(i)-x(j));
end;

h=h+l;
end
simplify(h);if(nargin == 3)
f = subs (h,'t',x0); %计算插值点的函数值
else
f=collect(h);
f = vpa(f,6); %将插值多项式的系数化成6位精度的小数
end结果:
f = 0.0201>>

⑧ matlab 拉格朗日插值法

function f = Language(x,y,x0)
%求已知数据点的拉格朗日插值多项式
%已知数据点的x坐标向量: x
%已知数据点的y坐标向量: y
%插值点的x坐标: x0
%求得的拉格朗日插值多项式或在x0处的插值: f

syms t;
if(length(x) == length(y))
n = length(x);
else
disp('x和y的维数不相等!');
return; %检错
end
f=0.0; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%改为 f=sym(0);
for(i=1:n)
l=y(i); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%改为 l=sym(y(i));
for(j=1:i-1)
l=l*(t-x(j))/(x(i)-x(j));
end;
for(j=i+1:n)
l=l*(t-x(j))/(x(i)-x(j));
end;

f=f+1;
simplify(f);

if(i==n)
if(nargin == 3)
f = subs (f,'t',x0); %计算插值点的函数值
else
f=collect(f);
f = vpa(f,6); %将插值多项式的系数化成6位精度的小数
end
end
end

-------------------------------------------------------
下面的这个应该可以:

function f = Language(x,y,x0)
%求已知数据点的拉格朗日插值多项式
%已知数据点的x坐标向量: x
%已知数据点的y坐标向量: y
%插值点的x坐标: x0
%求得的拉格朗日插值多项式或在x0处的插值: f

syms t l;
if(length(x) == length(y))
n = length(x);
else
disp('x和y的维数不相等!');
return; %检错
end

h=sym(0);
for (i=1:n)
l=sym(y(i));
for(j=1:i-1)
l=l*(t-x(j))/(x(i)-x(j));
end;
for(j=i+1:n)
l=l*(t-x(j))/(x(i)-x(j));
end;

h=h+l;
end
simplify(h);

if(nargin == 3)
f = subs (h,'t',x0); %计算插值点的函数值
else
f=collect(h);
f = vpa(f,6); %将插值多项式的系数化成6位精度的小数
end

阅读全文

与matlab拉格朗日插值算法相关的资料

热点内容
卡尔曼滤波算法书籍 浏览:763
安卓手机怎么用爱思助手传文件进苹果手机上 浏览:841
安卓怎么下载60秒生存 浏览:800
外向式文件夹 浏览:232
dospdf 浏览:428
怎么修改腾讯云服务器ip 浏览:382
pdftoeps 浏览:490
为什么鸿蒙那么像安卓 浏览:732
安卓手机怎么拍自媒体视频 浏览:183
单片机各个中断的初始化 浏览:721
python怎么集合元素 浏览:477
python逐条解读 浏览:829
基于单片机的湿度控制 浏览:496
ios如何使用安卓的帐号 浏览:879
程序员公园采访 浏览:807
程序员实战教程要多长时间 浏览:970
企业数据加密技巧 浏览:132
租云服务器开发 浏览:809
程序员告白妈妈不同意 浏览:332
攻城掠地怎么查看服务器 浏览:597