matlab高手请进,
来源:学生作业帮 编辑:作业帮 分类:综合作业 时间:2024/09/02 00:44:12
matlab高手请进,
function p=lag(n)
syms t;
p=0;
for k=1:(n+1)
x(k)=2*(k-1)/n-1;
end
for i=1:(n+1)
p1=fun(x(i));
for j=1:(n+1)
if j==i
p1=p1*1;
else
p1=p1*(t-x(j))/(x(i)-x(j));
end
end
p=p+p1;
end
p=collect(p);
function e=figerr(n)
format long
p=lag(n);
for j=1:5*n+1
x(j)=2*(j-1)/(5*n)-1;
y(j)=fun(x(j));
y1(j)=subs(p,'t',x(j));
e(j)=abs(y(j)-y1(j));
end
figure(1);
plot(x,y,'r');
hold on
plot(x,y1,'b');
grid on
figure(2);
s=1:5*n+1;
plot(s,e,'*r');
grid on
上面这个可以得出正确的图形
但用切比雪夫结点时,会出现问题,如下面的程序
function p=lag1(n)
syms t;
p=0;
for k=1:(n+1)
x(k)=cos((pi/2)*(2*k-1)/(n+1));
end
for i=1:(n+1)
p1=fun(x(i));
for j=1:(n+1)
if j==i
p1=p1*1;
else
p1=p1*(t-x(j))/(x(i)-x(j));
end
end
p=p+p1;
end
p=collect(p);
function e=figerr1(n)
format long
p=lag1(n);
for j=1:5*n+1
x(j)=2*(j-1)/(5*n)-1;
y(j)=fun(x(j));
y1(j)=subs(p,'t',x(j));
e(j)=abs(y(j)-y1(j));
end
figure(1);
plot(x,y,'r');
hold on
plot(x,y1,'b');
grid on
figure(2);
s=1:5*n+1;
plot(s,e,'*r');
grid on
![](http://img.wesiedu.com/upload/d/32/d32fada241ee2b10072407700374fe6e.jpg)
function p=lag(n)
syms t;
p=0;
for k=1:(n+1)
x(k)=2*(k-1)/n-1;
end
for i=1:(n+1)
p1=fun(x(i));
for j=1:(n+1)
if j==i
p1=p1*1;
else
p1=p1*(t-x(j))/(x(i)-x(j));
end
end
p=p+p1;
end
p=collect(p);
function e=figerr(n)
format long
p=lag(n);
for j=1:5*n+1
x(j)=2*(j-1)/(5*n)-1;
y(j)=fun(x(j));
y1(j)=subs(p,'t',x(j));
e(j)=abs(y(j)-y1(j));
end
figure(1);
plot(x,y,'r');
hold on
plot(x,y1,'b');
grid on
figure(2);
s=1:5*n+1;
plot(s,e,'*r');
grid on
上面这个可以得出正确的图形
但用切比雪夫结点时,会出现问题,如下面的程序
function p=lag1(n)
syms t;
p=0;
for k=1:(n+1)
x(k)=cos((pi/2)*(2*k-1)/(n+1));
end
for i=1:(n+1)
p1=fun(x(i));
for j=1:(n+1)
if j==i
p1=p1*1;
else
p1=p1*(t-x(j))/(x(i)-x(j));
end
end
p=p+p1;
end
p=collect(p);
function e=figerr1(n)
format long
p=lag1(n);
for j=1:5*n+1
x(j)=2*(j-1)/(5*n)-1;
y(j)=fun(x(j));
y1(j)=subs(p,'t',x(j));
e(j)=abs(y(j)-y1(j));
end
figure(1);
plot(x,y,'r');
hold on
plot(x,y1,'b');
grid on
figure(2);
s=1:5*n+1;
plot(s,e,'*r');
grid on
![](http://img.wesiedu.com/upload/d/32/d32fada241ee2b10072407700374fe6e.jpg)
我可是冲着那个500分来的
%by dynamic
%see also
%contact me matlabsky@gmail.com
%2009.2.21
%
function [c,L,str,yi]=lagrange(x,y,xi)
%Lagrange(拉格朗日)插值
%基于N+1个点,通过P(x)=sum(y(k)*L(N,k)(x),k=0...N) 其中k=0,1...,N 计算Lagrange多项式
%
%输入参数
%-x 横坐标数据
%-y 纵坐标数据
%-xi 需要插值的x坐标
%输出参数
%-c 拉格朗日插值多项式系数
%-L 拉格朗日系数多项式
%-str 拉格朗日差值多项式表达式
%-yi 插值结果
%Example
% x=[0, 1, 2, 3];
% y=[1, 3, 5, 12];
% xi=0:0.1:3;
% [c,L,str,yi]=lagrange(x,y,xi)
%
%by dynamic
%see also
%contact me matlabsky@gmail.com
%2009.2.21
%
w=length(x);
n=w-1;
L=zeros(w);
for kk=1:w
v=1;
for jj=1:w
if kk~=jj
v=conv(v,poly(x(jj)))/(x(kk)-x(jj));
end
end
L(kk,:)=v;
end
c=y*L;
str=poly2sym(c);
yi=polyval(x,xi);
plot(xi,yi,'r-.',x,y,'ob')
对于你的问题只要如下运行程序,就可以了
x=[0, 1, 2, 3];
y=[1, 3, 5, 12];
xi=0:0.1:3;
[c,L,str,yi]=lagrange(x,y,xi)
得到
c =
0.8333 -2.5000 3.6667 1.0000
L =
-0.1667 1.0000 -1.8333 1.0000
0.5000 -2.5000 3.0000 0
-0.5000 2.0000 -1.5000 0
0.1667 -0.5000 0.3333 0
str =
5/6*x^3-5/2*x^2+11/3*x+1
%by dynamic
%see also
%contact me matlabsky@gmail.com
%2009.2.21
%
function [c,L,str,yi]=lagrange(x,y,xi)
%Lagrange(拉格朗日)插值
%基于N+1个点,通过P(x)=sum(y(k)*L(N,k)(x),k=0...N) 其中k=0,1...,N 计算Lagrange多项式
%
%输入参数
%-x 横坐标数据
%-y 纵坐标数据
%-xi 需要插值的x坐标
%输出参数
%-c 拉格朗日插值多项式系数
%-L 拉格朗日系数多项式
%-str 拉格朗日差值多项式表达式
%-yi 插值结果
%Example
% x=[0, 1, 2, 3];
% y=[1, 3, 5, 12];
% xi=0:0.1:3;
% [c,L,str,yi]=lagrange(x,y,xi)
%
%by dynamic
%see also
%contact me matlabsky@gmail.com
%2009.2.21
%
w=length(x);
n=w-1;
L=zeros(w);
for kk=1:w
v=1;
for jj=1:w
if kk~=jj
v=conv(v,poly(x(jj)))/(x(kk)-x(jj));
end
end
L(kk,:)=v;
end
c=y*L;
str=poly2sym(c);
yi=polyval(x,xi);
plot(xi,yi,'r-.',x,y,'ob')
对于你的问题只要如下运行程序,就可以了
x=[0, 1, 2, 3];
y=[1, 3, 5, 12];
xi=0:0.1:3;
[c,L,str,yi]=lagrange(x,y,xi)
得到
c =
0.8333 -2.5000 3.6667 1.0000
L =
-0.1667 1.0000 -1.8333 1.0000
0.5000 -2.5000 3.0000 0
-0.5000 2.0000 -1.5000 0
0.1667 -0.5000 0.3333 0
str =
5/6*x^3-5/2*x^2+11/3*x+1