当前位置:文档之家› 第六章 函数的插值方法

第六章 函数的插值方法

第六章 函数的插值方法
第六章 函数的插值方法

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 求三个一次多项式)(x g 、)(x h 和)(x f 的积)()()(x h x g x f ??.它们的零点分别依次为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*eye(size (A)) - A)

6.2 拉格朗日(Lagrange )插值及其MATLAB 程序

6.2.1 线性插值及其MATLAB 程序

例 6.2.1 已知函数)(x f 在]3,1[上具有二阶连续导数,5)("≤x f ,且满足条件2)3(,1)1(==f f .求线性插值多项式和函数值)5.1(f ,并估计其误差.

解 输入程序

>> 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)

运行后输出基函数l 0和l 1及其插值多项式的系数向量P (略)、插值多项式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)))/2

运行后输出误差限为

R1 =

1.8750

例6.2.2 求函数=)(x f e x

-在]1,0[上线性插值多项式,并估计其误差.

解 输入程序

>> X=[0,1]; Y =exp(-X) ,

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),

运行后输出基函数l 0和l 1及其插值多项式的系数向量P 和插值多项式L 为

l0 = l1 = P =

-x+1 x -0.6321 1.0000

L =

-1423408956596761/2251799813685248*x+1

输入程序

>> M=1;x=0:0.001:1; R1=M*max(abs((x-X(1)).*(x-X(2))))./2

运行后输出误差限为

R1 =

0.1250.

6.2.2 抛物线插值及其MATLAB 程序

例 6.2.3 求将区间 [0, π/2] 分成n 等份)2,1(=n ,用x x f y cos )(==产生1+n 个节点,然后根据(6.9)和(6.13)式分别作线性插值函数)(1x P 和抛物线插值函数)(2x P .用它们分别计算cos (π/6) (取四位有效数字),并估计其误差.

解 输入程序

>> X=[0,pi/2]; Y =cos(X) ,

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=pi/6;

Y = polyval(P,x)

运行后输出基函数l 0和l 1及其插值多项式的系数向量P 、插值多项式和插值为

l0 =

-5734161139222659/9007199254740992*x+1

l1 =

5734161139222659/9007199254740992*x

P =

-0.6366 1.0000

L =

-5734161139222659/9007199254740992*x+1

Y =

0.6667

输入程序

>> M=1;x=pi/6; R1=M*abs((x-X(1))*(x-X(2)))/2

运行后输出误差限为

R1 =

0.2742.

(2) 输入程序

>> X=0:pi/4:pi/2; Y =cos(X) ,

poly(X(3)))/(( X(1)- X(2))* ( X(1)- X(3))),

l11= conv (poly(X(1)),

poly(X(3)))/(( X(2)- X(1))* ( X(2)- X(3))),

l21= conv (poly(X(1)),

poly(X(2)))/(( X(3)- X(1))* ( X(3)- X(2))),

l0=poly2sym (l01),l1=poly2sym (l11),l2=poly2sym (l21),

P = l01* Y(1)+ l11* Y(2) + l21* Y(3), L=poly2sym (P),x=pi/6;

Y = polyval(P,x)

运行后输出基函数l 01、l 11和l 21及其插值多项式的系数向量P 、插值多项式L 和插值Y 为

l0 =

228155022448185/281474976710656*x^2-2150310427208497/1125899906842624*x+1

l1 =

-228155022448185/140737488355328*x^2+5734161139222659/2251799813685248*x

l2 =

228155022448185/281474976710656*x^2-5734161139222659/9007199254740992*x

P =

-0.3357 -0.1092 1.0000

L=

-6048313895780875/18014398509481984*x^2-7870612110600739/72057594037927936

*x+1

Y =

0.8508

输入程序

>> M=1;x=pi/6; R2=M*abs((x-X(1))*(x-X(2)) *(x-X(3)))/6

运行后输出误差限为

R2 =

0.0239.

6.2.3 n 次拉格朗日(Lagrange )插值及其MATLAB 程序

例 6.2.4 给出节点数据00.17)00.2(=-f ,00.1)00.0(=f ,00.2)00.1(=f ,00.17)00.2(=f ,作三次拉格朗日插值多项式计算)6.0(f ,并估计其误差.

解 输入程序

>> X=[-2,0,1,2]; Y =[17,1,2,17];

p1=poly(X(1)); p2=poly(X(2));

p3=poly(X(3)); p4=poly(X(4));

l01= conv ( conv (p2, p3), p4)/(( X(1)- X(2))* ( X(1)- X(3))

* ( X(1)- X(4))),

l11= conv ( conv (p1, p3), p4)/(( X(2)- X(1))* ( X(2)- X(3))

* ( X(2)- X(4))),

l21= conv ( conv (p1, p2), p4)/(( X(3)- X(1))* ( X(3)- X(2))

* ( X(3)- X(4))),

l31= conv ( conv (p1, p2), p3)/(( X(4)- X(1))* ( X(4)- X(2))

* ( X(4)- X(3))),

l0=poly2sym (l01),

l1=poly2sym (l11),l2=poly2sym (l21), l3=poly2sym (l31),

P = l01* Y(1)+ l11* Y(2) + l21* Y(3) + l31* Y(4),

运行后输出基函数l 0,l 1,l 2和l 3及其插值多项式的系数向量P (略)为

l0 =

-1/24*x^3+1/8*x^2-1/12*x ,l1 =1/4*x^3-1/4*x^2-x+1

l2 =

-1/3*x^3+4/3*x ,l3 =1/8*x^3+1/8*x^2-1/4*x

输入程序

>> L=poly2sym (P),x=0.6; Y = polyval(P,x)

运行后输出插值多项式和插值为

L = Y =

x^3+4*x^2-4*x+1 0.2560.

输入程序

>> syms M; x=0.6;

R3=M*abs((x-X(1))*(x-X(2)) *(x-X(3)) *(x-X(4)))/24

运行后输出误差限为

R3 =

91/2500*M

R 3 )(04.0)4(ξ≈f , )00.2,00.2(-∈ξ.

6.2.5 拉格朗日多项式和基函数的MATLAB 程序

求拉格朗日插值多项式和基函数的MATLAB 主程序

function [C, L,L1,l]=lagran1(X,Y)

m=length(X); L=ones(m,m);

for k=1: m

V=1;

for i=1:m

if k~=i

V=conv(V,poly(X(i)))/(X(k)-X(i));

end

end

L1(k,:)=V; l(k,:)=poly2sym (V)

end

C=Y*L1;L=Y*l

例6.2.5 给出节点数据03.17)15.2(=-f ,24.7)00.1(=-f ,05.1)01.0(=f ,

03.2)02.1(=f ,

06.17)03.2(=f ,05.23)25.3(=f ,作五次拉格朗日插值多项式和基函数,并写出估计其误差的公式.

解 在MATLAB 工作窗口输入程序

>> 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];

[C, L ,L1,l]= lagran1(X,Y)

运行后输出五次拉格朗日插值多项式L 及其系数向量C ,基函数l 及其系数矩阵L 1如下

C =

-0.2169 0.0648 2.1076 3.3960 -4.5745 1.0954

L =

1.0954-4.5745*x+3.3960*x^2+

2.1076*x^3+0.0648*x^4-0.2169*x^5

L1 =

-0.0056 0.0299 -0.0323 -0.0292 0.0382 -0.0004

0.0331 -0.1377 -0.0503 0.6305 -0.4852 0.0048

-0.0693 0.2184 0.3961 -1.2116 -0.3166 1.0033

0.0687 -0.1469 -0.5398 0.6528 0.9673 -0.0097

-0.0317 0.0358 0.2530 -0.0426 -0.2257 0.0023

0.0049 0.0004 -0.0266 0.0001 0.0220 -0.0002

l =

[ -0.0056*x^5+0.0299*x^4-0.0323*x^3-0.0292*x^2+0.0382*x-0.0004]

[ 0.0331*x^5-0.1377*x^4-0.0503*x^3+0.6305*x^2-0.4852*x+0.0048]

[ -0.0693*x^5+0.2184*x^4+0.3961*x^3-1.2116*x^2-0.3166*x+1.0033]

[ 0.0687*x^5-0.1469*x^4-0.5398*x^3+0.6528*x^2+0.9673*x-0.0097]

[ -0.0317*x^5+0.0358*x^4+0.2530*x^3-0.0426*x^2-0.2257*x+0.0023]

[ 0.0049*x^5+0.0004 *x^4-0.0266*x^3+0.0001*x^2+0.0220*x-0.0002]

估计其误差的公式为

)(5x R )25.3)(03.2)(02.1)(01.0()00.1)(15.2(!

6)()6(----++=x x x x x x f ξ,)3.25,-2.15(∈ξ.

6.2.6 拉格朗日插值及其误差估计的MATLAB 程序

拉格朗日插值及其误差估计的MATLAB 主程序

function [y,R]=lagranzi(X,Y,x,M)

n=length(X); m=length(x);

for i=1:m

z=x(i);s=0.0;

for k=1:n

p=1.0; q1=1.0; c1=1.0;

for j=1:n

if j~=k

p=p*(z-X(j))/(X(k)-X(j));

end

q1=abs(q1*(z-X(j)));c1=c1*j;

end

s=p*Y(k)+s;

end

y(i)=s;

end

R=M*q1/c1;

例 6.2.6 已知5.030sin = ,1707.045sin = ,0866.060sin = ,用拉格朗

日插值及其误差估计的MATLAB 主程序求

40sin 的近似值,并估计其误差.

解 在MATLAB 工作窗口输入程序

>> x=2*pi/9; M=1; X=[pi/6 ,pi/4, pi/3];

Y=[0.5,0.7071,0.8660]; [y,R]=lagranzi(X,Y,x,M)

运行后输出插值y 及其误差限R 为

y = R =

0.6434 8.8610e-004.

6.3 牛顿(Newton )插值及其MATLAB 程序

6.3.3 牛顿插值多项式、差商和误差公式的MATLAB 程序

求牛顿插值多项式和差商的MATLAB 主程序

function [A,C,L,wcgs,Cw]= newploy(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:n

for i=j:n

A(i,j)=(A(i,j-1)- A(i-1,j-1))/(X(i)-X(i-j+1));

end

b=poly(X(j-1));q1=conv(q,b); c1=c1*j; q=q1;

end

C=A(n,n); b=poly(X(n)); q1=conv(q1,b);

for k=(n-1):-1:1

C=conv(C,poly(X(k))); d=length(C); C(d)=C(d)+A(k,k);

end

L(k,:)=poly2sym(C); Q=poly2sym(q1);

syms M

wcgs=M*Q/c1; Cw=q1/c1;

例6.3.3 给出节点数据03.17)15.2(=-f ,24.7)00.1(=-f ,05.1)01.0(=f ,03.2)02.1(=f ,06.17)03.2(=f ,05.23)25.3(=f ,作五阶牛顿插值多项式和差商,并写出其估计误差的公式.

解 (1)保存名为newpoly.m 的M 文件.

(2)输入MATLAB 程序

Y=[17.03 7.24 1.05 2.03 17.06 23.05];

[A,C,L,wcgs,Cw]= newdcw (X,Y)

运行后输出差商矩阵A ,五阶牛顿插值多项式L 及其系数向量C , 插值余项公式L 及其向量C w 如下

A =

17.0300 0 0 0 0 0

7.2400 -8.5130 0 0 0 0

1.0500 -6.1287 1.1039 0 0 0

2.0300 0.9703

3.5144 0.7604 0 0

17.0600 14.8812 6.8866 1.1129 0.0843 0

23.0500 4.9098 -4.4715 -3.5056 -1.0867 -0.2169

C =

-0.2169 0.0648 2.1076 3.3960 -4.5745 1.0954

L =

-7813219284746629/36028797018963968*x^5+583849564517807/9007199

254740992*x^4+593245028711263/281474976710656*x^3+3823593773002357/1125899906842624*x^2-321902673270315/70368744177664*x+308328649211299/281474976710656

wcgs =

1/720*M*(x^6-79/25*x^5-14201/2500*x^4+4934097026900981/2814749767

10656*x^3+154500712237335/35184372088832*x^2-8170642380559269/562949953421312*x+5212760744134241/36028797018963968)

Cw =

0.0014 -0.0044 -0.0079 0.0243 0.0061 -0.0202 0.0002

L =1.0954-4.5745*x +3.3960*x ^2+2.1076*x ^3+0.0648*x ^4-0.2169*x ^5.

估计其误差的公式为

)(5x R )25.3)(03.2)(02.1)(01.0()00.1)(15.2(!

6)()6(----++=x x x x x x f ξ,)3.25,-2.15(∈ξ.

例6.3.4 求函数7)(-=x f e 5/x -在]6,2[上五阶牛顿插值多项式,估计其误差的公式和误差限公式.用它们计算)156.3(f ,并估计其误差.

解 (1)输入MATLAB 程序

>> X=2:4/5:6; Y=-7*exp(-X/5);

[A,C,L,wcgs,Cw]= newploy(X,Y), x1=2:0.001:6;

M=max(-7*exp(-x1/5)/(5^6)),

运行后输出差商矩阵A , 五阶牛顿插值多项式L 及其系数向量C , 插值余项公式L 及其向量C w 如下

A =

-4.6922 0 0 0 0 0

-3.9985 0.8672 0 0 0 0

-3.4073 0.7390 -0.0801 0 0 0

-2.9035 0.6297 -0.0683 0.0049 0 0

-2.4742 0.5366 -0.0582 0.0042 -0.0002 0

-2.1084 0.4573 -0.0496 0.0036 -0.0002 0.0000

C =

0.0000 -0.0004 0.0089 -0.1389 1.3985 -6.9991

L =

9721799720875/1152921504606846976*x^5-3503994098647815/9223

372036854775808*x^4+160742008798419/18014398509481984*x^3-1251152213853501/9007199254740992*x^2+6298131904328647/4503599627370496*x-3940156929554013/562949953421312

wcgs =

1/720*M*(x^6-24*x^5+1172/5*x^4-5952/5*x^3+7276634802928539/219

9023255552*x^2-5237461186650519/1099511627776*x+6085939356447121/219902

3255552)

Cw =

0.0014 -0.0333 0.3256 -1.6533 4.5959 -6.6159 3.8438

M =

-1.3494e-004

(2)输入MATLAB 程序

>> syms x

wcgs1=1/720*M*1/720*M*(x^6-24*x^5+1172/5*x^4-5952/5*x^3

+7276634802928539/2199023255552*x^2-5237461186650519/1099511627776*x+6085939356447121/2199023255552),

运行后输出误差限公式wcgs1如下

wcgs1 =

5565367633581443/158456325028528675187087900672*x^6-16696102900744329/19807040628566084398385987584*x^5+1630652716639362799/198070406285660843983859875840*x^4-517579189923074199/12379400392853802748991242240*x^3+40497147813610772932297365501777/348449143727040986586495598010130648530944*x^2-29148396970323855270001164718917/174224571863520493293247799005065324265472*x+33870489914310283926665276375603/348449143727040986586495598010130648530944

(3)输入MATLAB 程序

>> x=3.156;

y=9721799720875/1152921504606846976*x^5-3503994098647815

/9223372036854775808*x^4+160742008798419/18014398509481984*x^3-1251152213853501/9007199254740992*x^2+6298131904328647/4503599627370496*x-3940156929554013/562949953421312

wcgs2=1/720*M*(x^6-24*x^5+1172/5*x^4-5952/5*x^3+72766348

02928539/2199023255552*x^2-5237461186650519/1099511627776*x+6085939356447121/2199023255552)

运行后输出)156.3(f 的近似值y ,及其误差限wcgs 2如下

y =

-3.7237

wcgs2 =

-2.4764e-007

6.3.4 牛顿插值及其误差估计的MATLAB 程序

牛顿插值及其误差估计的MATLAB 主程序

function [y,R]= newcz(X,Y,x,M)

n=length(X); m=length(x);

for t=1:m

z=x(t); A=zeros(n,n);A(:,1)=Y';

s=0.0; p=1.0; q1=1.0; c1=1.0;

for j=2:n

for i=j:n

A(i,j)=(A(i,j-1)- A(i-1,j-1))/(X(i)-X(i-j+1));

end

q1=abs(q1*(z-X(j-1)));c1=c1*j;

end

C=A(n,n);q1=abs(q1*(z-X(n)));

for k=(n-1):-1:1

C=conv(C,poly(X(k)));d=length(C); C(d)=C(d)+A(k,k);

end

y(k)= polyval(C, z);

end

R=M*q1/c1;

例6.3.5 已知5.030sin =

,1707.045sin = ,0866.060sin = ,用牛顿插值

法求

40sin 的近似值,估计其误差,并与例 6.2.6的计算结果比较.

解 方法1(牛顿插值及其误差估计的MATLAB 主程序)输入MATLAB 程序

>> x=2*pi/9;M=1; X=[pi/6 ,pi/4, pi/3];

Y=[0.5,0.7071,0.8660]; [y,R]= newcz(X,Y,x,M)

运行后输出

y = R =

0.6434 8.8610e-004

方法2 (求牛顿插值多项式和差商的MATLAB 主程序)输入MATLAB 程序

>> x=2*pi/9; X=[pi/6 ,pi/4, pi/3];

Y=[0.5,0.7071,0.8660]; M=1;

[A,C,L,wcgs,Cw]= newploy(X,Y), y=polyval(C,x)

运行后输出结果

A =

0.5000 0 0

0.7071 0.7911 0

0.8660 0.6070 -0.3516

C =

-0.3516 1.2513 -0.0588

L =

-1583578379808357/4503599627370496*x^2+1408883391907005/

1125899906842624*x-132405829044691/2251799813685248

wcgs =

1/6*M*(x^3-3/4*x^2*pi+4012734077357799/2251799813685248*

x-7757769783530263/18014398509481984)

Cw =

0.1667 -0.3927 0.2970 -0.0718

y =

0.6434

上述两种方法计算y 的结果相同.

6.3.5 牛顿插值法的MATLAB 综合程序

求牛顿插值多项式、差商、插值及其误差估计的MATLAB 主程序

function [y,R,A,C,L]=newdscg(X,Y,x,M)

n=length(X); m=length(x);

for t=1:m

z=x(t); A=zeros(n,n);A(:,1)=Y';

s=0.0; p=1.0; q1=1.0; c1=1.0;

for j=2:n

for i=j:n

A(i,j)=(A(i,j-1)- A(i-1,j-1))/(X(i)-X(i-j+1));

end

q1=abs(q1*(z-X(j-1)));c1=c1*j;

end

C=A(n,n);q1=abs(q1*(z-X(n)));

for k=(n-1):-1:1

C=conv(C,poly(X(k)));

d=length(C);C(d)=C(d)+A(k,k);

end

y(k)= polyval(C, z);

end

R=M*q1/c1;L(k,:)=poly2sym(C);

例6.3.6 给出节点数据00.27)00.4(=-f ,00.1)00.0(=f ,00.2)00.1(=f ,00.17)00.2(=f ,作三阶牛顿插值多项式,计算)345.2(-f ,并估计其误差.

解 首先将名为newdscg.m 的程序保存为M 文件,然后在MA TLAB 工作窗口输入程序

>> syms M,X=[-4,0,1,2]; Y =[27,1,2,17]; x=-2.345;

[y,R,A,C,P]=newdscg(X,Y,x,M)

运行后输出插值y )345.2(-≈f 及其误差限公式R ,三阶牛顿插值多项式P 及其系数向量C ,差商的矩阵A 如下

y =

22.3211

R =

1323077530165133/562949953421312*M (即R =2.3503*M )

A=

27.0000 0 0 0

1.0000 -6.5000 0 0

2.0000 1.0000 1.5000 0

17.0000 15.0000 7.0000 0.9167

C =

0.9167 4.2500 -4.1667 1.0000

P =

11/12*x^3+17/4*x^2-25/6*x+1

例 6.3.7 将区间 [0,π/2] 分成n 等份)3,2(=n ,用x x f y sin )(==产生1+n 个节点,求二阶和三阶牛顿插值多项式)(2x P 和)(3x P .用它们分别计算sin (π/7) (取四位有效数字),并估计其误差.

解 首先将名为newdscg.m 的程序保存为M 文件,然后在MA TLAB 工作窗口输入程序

>> X1=0:pi/4:pi/2; Y1 =sin(X1); M=1; x=pi/7;

X2=0:pi/6:pi/2; Y2 =sin(X2);

[y1,R1,A1,C1,P2]=newdscg(X1,Y1,x,M),

[y2,R2,A2,C2,P3]=newdscg(X2,Y2,x,M)

运行后输出插值y 1=)7/sin()7/(2π≈πP 和y 2=)7/sin()7/(3π≈πP 及其误差限R 1和R 2,二阶和三阶牛顿插值多项式P 2和P 3及其系数向量C 1和C 2,差商的矩阵A 1和A 2如下

y1 =

0.4548

R1 =

0.0282

A1 =

0 0 0

0.7071 0.9003 0

1.0000 0.3729 -0.3357

C1 =

-0.3357 1.1640 0

P2 =

-3024156947890437/9007199254740992*x^2+163820246322191/140737488355328*x

y2 =

0.4345

R2 =

9.3913e-004

A2 =

0 0 0 0

0.5000 0.9549 0 0

0.8660 0.6991 -0.2443 0

1.0000 0.2559 -0.4232 -0.1139

C2 =

-0.1139 -0.0655 1.0204 0

P3 =

-1025666884451963/9007199254740992*x^3-4717668559161127/7

2057594037927936*x^2+4595602396951547/4503599627370496*x

6.4 埃尔米特(Hermite )插值及其MATLAB 程序

6.4.3 埃尔米特插值多项式和误差公式的MATLAB 程序

求埃尔米特插值多项式和误差公式的MATLAB 主程序

function [Hc, Hk,wcgs,Cw]= hermite (X,Y,Y1)

m=length(X); n=M1;s=0; H=0;q=1;c1=1; L=ones(m,m);

G=ones(1,m);

for k=1:n+1

V=1;

for i=1:n+1

if k~=i

s=s+(1/(X(k)-X(i)));

V=conv(V,poly(X(i)))/(X(k)-X(i));

end

h=poly(X(k)); g=(1-2*h*s); G=g*Y(k)+h*Y1(k);

end

H=H+conv(G,conv(V,V)); b=poly(X(k));b2=conv(b,b);

q=conv(q,b2); t=2*n+2;

Hc=H;Hk=poly2sym (H); Q=poly2sym(q);

end

for i=1:t

c1=c1*i;

end

syms M,wcgs=M*Q/c1; Cw=q/c1;

例6.4.3 给定函数)(x f 在点2/,4/,6/210π=π=π=x x x 处的函数值5.0)(0=x f , 1707.0)(1=x f ,1)(2=x f 和导数值0866.0)(0'=x f ,1707.0)(1'=x f ,

0000.0)(0'=x f ,且1)()6(≤x f ,求函数)(x f 在点210,,x x x 处的五阶埃尔米特插值

多项式)(5x H 和误差公式,计算)1.567

(f 并估计其误差. 解 (1)保存名为hermite.m 的M 文件.

(2)在MATLAB 工作窗口输入程序

>>X=[pi/6,pi/4,pi/2]; Y=[0.5,0.7071,1];

Y1=[0.8660,0.7071,0]; [Hc, Hk,wcgs,Cw]= hermite (X,Y,Y1)

运行后输出五阶埃尔米特插值多项式H k 及其系数向量H c ,误差公式wcgs 及其系数向量C w 如下

Hc =

1.0e+003 *

0.1912 -0.9273 1.6903 -1.4380 0.5751 -0.0866

Hk =

6725828781679091/35184372088832*x^5-4078286086775209/4398

046511104*x^4+7434035571017927/4398046511104*x^3-3162205449085973/2199023255552*x^2+5058863928652835/8796093022208*x-6094057839958843/70368744177664

wcgs =

1/720*M*(x^6-11/6*x^5*pi+7446708432019761/562949953421312

*x^4-4363745503235773/281474976710656*x^3+21569239021155/2199023255552*x^2-7178073637328281/2251799813685248*x+3758430567659515/9007199254740992)

Cw =

0.0014 -0.0080 0.0184 -0.0215 0.0136 -0.0044 0.0006 当M x f =≤1)()6(时的误差公式为

R=0.001 4* x ^6-0.008 0*x ^5+0.018 4*x ^4-0.021 5*x ^3+0.013 6*x ^2-0.004 4*x +0.000 6

(3)在MATLAB 工作窗口输入程序

>> x=1.567;M=1;

Hk=6725828781679091/35184372088832*x^5-4078286086775209/

4398046511104*x^4+7434035571017927/4398046511104*x^3-3162205449085973/2199023255552*x^2+5058863928652835/8796093022208*x-6094057

wcgs=1/720*M*(x^6-11/6*x^5*pi+7446708432019761/562949953

421312*x^4-4363745503235773/281474976710656*x^3+21569239021155/2199023255552*x^2-7178073637328281/2251799813685248*x+3758430567659515/9007199254740992)

运行后输出)1.567

(f 的近似值Hk 及其误差wcgs 如下 Hk =

2.5265

wcgs =

1.3313e-008

6.5 分段插值及其MATLAB 程序

6.5.1 高次插值的振荡和MATLAB 程序

例6.5.1 作下列函数在指定区间上的n 次拉格朗日插值多项式)(x L n

)10,8,6,4,2(=n 的图形,并讨论插值多项式)(x L n 的次数与误差)(x R n 的关系.

(1)))432sin 3tan(cos()(2

x x x f ++=,],[ππ-∈x ;(2)211)(x x f +=,]5,5[-∈x . 解 将计算n 次拉格朗日插值多项式)(x L n 的值的MATLAB 程序保存名为

lagr1.m 的M 文件.

function y=lagr1(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

(1)现提供作n 次拉格朗日插值多项式)(x L n 的图形的MATLAB 程序,保存名为

Li1.m 的M 文件

m=150; x=-pi:2*pi/(m-1): pi;

y=tan(cos((3^(1/2)+sin(2*x))./(3+4*x.^2)));

plot(x,y,'k-'),

gtext('y=tan(cos((sqrt(3)+sin(2x))/(3+4x^2)))'),pause

n=3; x0=-pi:2*pi/(3-1):pi;

y0=tan(cos((3^(1/2)+sin(2*x0))./(3+4*x0.^2)));

y1=lagr1(x0,y0,x);hold on,

plot(x,y1,'g<'),gtext('n=2'),pause,hold off

n=5; x0=-pi:2*pi/(n-1):pi;

y0=tan(cos((3^(1/2)+sin(2*x0))./(3+4*x0.^2)));

y2=lagr1(x0,y0,x);hold on,

plot(x,y2,'b:'),gtext('n=4'),pause,hold off

n=7; x0=-pi:2*pi/(n-1):pi;

y0=tan(cos((3^(1/2)+sin(2*x0))./(3+4*x0.^2)));

y3=lagr1(x0,y0,x);hold on,

plot(x,y3,'rp'),gtext('n=6'),pause,hold off

n=9; x0=-pi:2*pi/(n-1):pi;

y4=lagr1(x0,y0,x);hold on,

plot(x,y4,'m*'),gtext('n=8'),pause,hold off

n=11; x0=-pi:2*pi/(n-1):pi;

y0=tan(cos((3^(1/2)+sin(2*x0))./(3+4*x0.^2)));

y5=lagr1(x0,y0,x);hold on,

plot(x,y5,'g:'),gtext('n=10')

title('高次拉格朗日插值的振荡')

在MA TLAB 工作窗口输入名为 Li1.m 的M 文件的文件名

>> Li1.m

回车运行后,便会逐次画出)(x f 在区间],[ππ-上的n 次拉格朗日插值多项式)(x L n )10,8,6,4,2(=n 的图形(略).

(2)在MATLAB 工作窗口输入程序

m=101; x=-5:10/(m-1):5; y=1./(1+x.^2);

z=0*x;plot(x,z,'r',x,y,'k-'),

gtext('y=1/(1+x^2)'),pause

n=3; x0=-5:10/(n-1):5; y0=1./(1+x0.^2);

y1=lagr1(x0,y0,x);hold on,

plot(x,y1,'g'),gtext('n=2'),pause,hold off

n=5; x0=-5:10/(n-1):5; y0=1./(1+x0.^2);

y2=lagr1(x0,y0,x);hold on,

plot(x,y2,'b:'),gtext('n=4'),pause,hold off

n=7; x0=-5:10/(n-1):5; y0=1./(1+x0.^2);

y3=lagr1(x0,y0,x);hold on,

plot(x,y3,'r'),gtext('n=6'),pause,hold off

n=9; x0=-5:10/(n-1):5; y0=1./(1+x0.^2);

y4=lagr1(x0,y0,x);hold on,

plot(x,y4,'r:'),gtext('n=8'),pause,hold off

n=11; x0=-5:10/(n-1):5; y0=1./(1+x0.^2);

y5=lagr1(x0,y0,x);hold on,

plot(x,y5,'m'),gtext('n=10')

title('高次拉格朗日插值的振荡')

回车运行后,便会逐次画出)1/(12x y +=在区间]5,5[-上的n 次拉格朗日插值多项式

)(x L n )10,8,6,4,2(=n 的图形(略).

6.5.4 作有关分段线性插值图形的MATLAB 程序

作有关分段线性插值图形的MATLAB 主程序

function s= xxczhjt1(x0,y0,xi,x,y)

s= interp1(x0,y0,xi);

Sn= interp1(x0,y0,x0);

plot(x0,y0,'o',x0,Sn,'-',xi,s,'*',x,y,'-.')

legend('节点(xi,yi)','分段线性插值函数Sn(x)','插值点(x,s)','

被插值函数y')

我们也可以直接在在MA TLAB 工作窗口编程序.例如,

>> x0 =-6:6; y0 =sin(x0);

xi = -6:.25:6;

yi = interp1(x0,y0,xi);

x=-6:0.001:6; y=sin(x);plot(x0,y0,'o',xi,yi,x,y,':'),

legend('节点(xi,yi)','分段线性插值函数','被插值函数y=sinx ')

title('y=sinx 及其分段线性插值函数和节点的图形')

>> x0 =-6:6; y0 =cos(x0);

xi = -6:.25:6;

yi = interp1(x0,y0,xi);

'o'':'legend('节点(xi,yi)','分段线性插值函数','被插值函数y=cosx ')

title('y=cosx 及其分段线性插值函数和节点的图形')

例6.5.5 设函数22511)(x

x f +=,在区间]1,1[-上取等距节点),(i i y x 10,,2,1,0 =i , 构造分段线性插值函数)(x S n ,用MA TLAB 程序计算各小区间的中点i x 处)(x S n 的值,作出节点,插值点,)(x f 和)(x S n 的图形.

解 节点的横坐标和插值点等取值与例6.5.4相同.在MATLAB 工作窗口输入程序

>>x0=-1:0.2:1; y0=1./(1+25.*x0.^2);

xi=-0.9:0.2:0.9;

b=max(x0);

a=min(x0);x=a:0.001:b;

y=1./(1+25.*x.^2);

s=xxczhjt1(x0,y0,xi,x,y), title('y=1/(1+25 x^2)的分段线性插

值的有关图形')

运行后屏幕显示各小区间中点i x 处)(x S n 的值,出现节点,插值点,)(x f 和)(x S n 的图形(略)

s =

Columns 1 through 4

0.04864253393665 0.07941176470588 0.15000000000000 0.35000000000000

Columns 5 through 8

0.75000000000000 0.75000000000000 0.35000000000000 0.15000000000000

Columns 9 through 10

0.07941176470588 0.04864253393665

6.5.5 用MATLAB 计算有关分段线性插值的误差

例6.5.8 设函数))432sin 3tan(cos()(2x

x x f ++=,在区间],[ππ-上取等距节点),(i i y x ,9,,2,1,0 =i ,构造分段线性插值函数)(x S n .

(1)用MA TLAB 程序计算各小区间中点i x 处)(x S n 的值,作出节点,插值点,)(x f 和)(x S n 的图形;

(2) 用MA TLAB 程序计算各小区间中点处)(x S n 的值及其相对误差;

(3) 用MA TLAB 程序估计)(max "

x f x ππ≤≤-和)(x S n 在区间],[ππ-上的误差限. 解 在MATLAB 工作窗口输入程序

>>h=2*pi/9; x0=-pi:h:pi;

y0=tan(cos((sqrt(3)+sin(2*x0))./(3+4*x0.^2)));

xi=-pi+h/2:h:pi-h/2;

fi=tan(cos((3^(1/2)+sin(2*xi))./(3+4*xi.^2)));

b=max(x0); a=min(x0);

x=a:0.001:b; y=tan(cos((sqrt(3)+sin(2.*x))./(3+4*x.^2))); si=xxczhjt1(x0,y0,xi,x,y);

Ri=abs((fi-si)./fi);

xi,fi,si,Ri,i=[xi',fi',si',Ri']

title('y=tan(cos((sqrt(3)+sin(2 x))/(3+4 x^2)))的分段线性插

值的有关图形')

运行后屏幕显示R i (略),并且作出节点,插值点,)(x f 和)(x S n 的图形(略).

在MATLAB 工作窗口输入程序

y=tan(cos((3^(1/2)+sin(2*x))/(3+4*x^2)));yxx=diff(y,x,2), 运行后屏幕显示函数)(x f 的二阶导数)(''x f (略).

在MATLAB 工作窗口输入程序

>> h=2*pi/9; x=-pi:0.0001:-pi;

yxx=2*tan(cos((3^(1/2)+sin(2.*x))./(3+4.*x.^2))).*(1+tan(

cos((3.^(1/2)+sin(2.*x))./(3+4.*x.^2))).^2).*sin((3^(1/2)+sin(2.*x))./(3+4.*x.^2)).^2.*(2*cos(2.*x)./(3+4*x.^2)-8*(3^(1/2)+sin

(2.*x))./(3+4.*x.^2).^2.*x).^2-(1+tan(cos((3^(1/2)+sin(2.*x))./(3+4.*x.^2))).^2).*cos((3^(1/2)+sin(2.*x))./(3+4.*x.^2)).*(2.*c os(2.*x)./(3+4.*x.^2)-8*(3^(1/2)+sin(2.*x))./(3+4.*x.^2).^2.*x).^2-(1+tan(cos((3^(1/2)+sin(2.*x))./(3+4.*x.^2))).^2).*sin((3.^(1/2)+sin(2.*x))./(3+4.*x.^2)).*(-4.*sin(2.*x)./(3+4.*x.^2)-32*cos(2.*x)./(3+4.*x.^2).^2.*x+128*(3^(1/2)+sin(2.*x))./(3+4.*x.^

2).^3.*x.^2-8*(3^(1/2)+sin(2.*x))./(3+4.*x.^2).^2);

myxx=max(yxx), R=(h^2)* myxx/8

运行后屏幕显示myxx =)(max "

x f x π≤≤π-和)(x S n 在区间],[ππ-上的误差限R 如下 myxx = R =

-0.02788637150664 -0.00169893490711

6.6 分段埃尔米特插值及其MATLAB 程序

6.6.2 分段埃尔米特插值的MATLAB 程序

调用格式一:YI=interp1(X,Y ,XI,'pchip')

调用格式二:YI=interp1(X,XI,'pchip')

例6.6.5 试用MA TLAB 程序计算例6.6.1中在各小区间中点处分段三次埃尔米特插值)(2/1+i n x H 及其相对误差.

解 在MATLAB 工作窗口输入程序

>> h=0.2;x0=-1:h:1;y0=1./(1+25.*x0.^2); xi=-0.9:h:0.9;

fi=1./(1+25.*xi.^2); yi=interp1(x0,y0,xi,'pchip');

Ri=abs((fi-yi)./fi); xi,fi,yi,Ri,i=[xi',fi',yi',Ri'] 运行后屏幕显示各小区间中点x i 处的函数值f i ,插值s i ,相对误差值R i (略).

6.6.3 作有关分段埃尔米特插值图形的MATLAB 程序

作有关分段埃尔米特插值图形的MATLAB 主程序

function H=hermitetx(x0,y0,xi,x,y)

H= interp1(x0,y0,xi,'pchip');

Hn= interp1(x0,y0,x,'pchip');

plot(x0,y0,'o',x,Hn,'-',xi,H,'*',x,y,'-.')

legend('节点(xi,yi)', '分段埃尔米特插值函数','插值点(x,H)','被

插值函数y')

我们也可以直接在在MATLAB 工作窗口编程序,例如,

>> x0 =-6:6; y0 =sin(x0); xi = -6:.25:6;

yi = interp1(x0,y0,xi,'pchip');

x=-6:0.001:6; y=sin(x); plot(x0,y0,'o',xi,yi,x,y,':'),

legend('节点(xi,yi)','分段埃尔米特插值函数','被插值函数

y=sinx')

title(' y=sinx 及其分段埃尔米特插值函数和节点的图形')

>> x0 =-6:6; y0 =cos(x0);

xi = -6:.25:6;yi = interp1(x0,y0,xi,'pchip');

x=-6:0.001:6; y=cos(x); plot(x0,y0,'o',xi,yi,x,y,':'),

legend('节点(xi,yi)','分段埃尔米特插值函数','被插值函数

y=cosx')

title(' y=cosx 及其分段埃尔米特插值函数和节点的图形')

例6.6.6 设函数2

11)(x x f +=定义在区间]5,5[-上,节点(X (i ),f (X (i )))的横坐标向量X 的元素是首项a =-5,末项b =5,公差h =1.5的等差数列,构造三次分段埃尔米特插值函数)(3,x H n .把区间]5,5[-分成20等份,构成20个小区间,用MA TLAB 程序计算各小区间中点i x 处)(3,x H n 的值,并作出节点,插值点,)(x f 和)(3,x H n 的图形.

解 在MATLAB 工作窗口输入程序

>>x0=-5:1.5:5;

y0=1./(1+x0.^2); x1=-4.75:0.5:4.75;

x=-5:0.001:5; y=1./(1+x.^2); H= hermitetx(x0,y0,x1,x,y)

title('函数y=1/(1+x^2)及其分段埃尔米特插值函数,插值,节点(xi,yi)

的图形')

运行后屏幕显示各小区间中点i x 处)(3,x H n 的值,出现节点,插值点,)(x f 和

)(3,x H n 的图形(略).

例6.6.7 设函数x x x f cos 5.0)(-=定义在区间],[ππ-上,取7=n ,按等距节点构造分段埃尔米特插值函数)(3,7x H ,用MATLAB 程序计算各小区间中点i x 处)(3,7x H 的值,作出节点,插值点,)(x f 和)(3,7x H 的图形.

解 记节点的横坐标7,,2,1,0,7/2, =π=+π-=i h ih x i 插值点)(2

1121

+++=i i i x x x ,6,,2,1,0 =i .在MATLAB 工作窗口输入程序 >> h=2*pi/7; x0=-pi:h:pi;

y0=0.5.*x0-cos(x0); xi=-pi+h/2:h:pi-h/2;

b=max(x0); a=min(x0); x=a:0.001:b;

y=0.5.*x-cos(x); H= hermitetx(x0,y0,xi,x,y)

title('函数y=0.5x-cos(x)及其分段埃尔米特插值函数,插值,节点

(xi,yi) 的图形')

运行后屏幕显示各小区间中点i x 处)(3,7x H 的值,出现节点,插值点,)(x f 和)(3,7x H 的图形(略).

6.6.4 用MATLAB 计算有关分段埃尔米特插值的误差

例6.6.8 设函数22511)(x

x f +=定义在区间]1,1[-上,取10=n ,按等距节点构造分段埃尔米特插值函数)(3,x H n ,用MATLAB 程序在]1,1[-上计算)(max )4(1

1x f x ≤≤-和)(3,x H n 的误差公式和误差限.

解 在MATLAB 工作窗口输入程序

>> syms h,x=-1:0.0001:1;

yxxxx=150000000./(1+25.*x.^2).^5.*x.^4-4500000./(1+25.*x.

^2).^4.*x.^2+15000./(1+25.*x.^2).^3;

myxxxx=max(yxxxx), R=(h^4)* abs(myxxxx/384)

运行后输出)(x f 的4阶导数在区间]1,1[-上绝对值的最大值myxxxx 和)(3,x H n 在区间

]1,1[-上的误差公式myxxxx 为

myxxxx = R =

15000 625/16*h^4

(4)在MATLAB 工作窗口输入程序

运行后输出误差限为

R =

0.06250000000000

例6.6.9 设函数))432sin 3tan(cos(

)(2x

x x f ++=定义在区间],[ππ-上,取9=n ,按等距节点构造分段埃尔米特插值函数)(3,x H n . (1)用MA TLAB 程序计算各小区间中点i x 处)(3,x H n 的值,作出节点,插值点,)(x f 和)(3,x H n 的图形;

(2) 并用MATLAB 程序计算各小区间中点处)(3,x H n 的值及其相对误差;

(3) 用MA TLAB 程序求)(max )4(x f x ππ≤≤-和)(3,x H n 在区间],[ππ-上的误差公式和各插值的误差限.

解 (1)记节点的横坐标9,,2,1,0,9/2, =π=+π-=i h ih x i ,插值点)(2

1121+++=i i i x x x ,8,,2,1,0 =i .在MATLAB 工作窗口输入程序

>>h=2*pi/9; x0=-pi:h:pi;

y0=tan(cos((sqrt(3)+sin(2*x0))./(3+4*x0.^2)));

xi=-pi+h/2:h:pi-h/2;

fi=tan(cos((3^(1/2)+sin(2*xi))./(3+4*xi.^2)));

b=max(x0); a=min(x0); x=a:0.001:b;

y=tan(cos((3^(1/2)+sin(2.*x))./(3+4*x.^2)));

Hi= hermite tx (x0,y0,xi,x,y);

Ri=abs((fi-Hi)./fi); xi,fi,Hi,Ri,i=[xi',fi',Hi',Ri']

title('函数y=tan(cos((sqrt(3)+sin(2x))/(3+4x^2)))及其分段埃

尔米特插值函数,插值,节点(xi,yi) 的图形')

运行后屏幕显示各小区间中点x i 处的函数值f i ,插值H i ,相对误差值R i ,并且作出节点,插值点,)(x f 和)(3,x H n 的图形(略).

(2)在MATLAB 工作窗口输入程序

>> syms x

y=tan(cos((3^(1/2)+sin(2*x))/(3+4*x^2)));

yxxxx=diff(y,x,4),%simple(yxxxx)

运行后屏幕显示函数)(x f 的4阶导数)()4(x f ,然后将输出的)()4(x f

编程求)(max )4(x f x ππ≤≤-和)(3,x H n 及其在区间],[ππ-上的误差限的MA TLAB 程序如下

>>syms h,x=-pi:0.0001:pi;

yxxxx=-12.*(1.+tan(cos((3.^(1./2)+sin(2.*x))./(3.+4.*x.^

2))).^2).^2.*sin((3.^(1./2)+sin(2.*x))./(3.+4.*x.^2)).^3.*(2.*c os(2.*x)./(3.+4.*x.^2)-8.*(3.^(1./2)+sin(2.*x))./(3.+4.*x.^2).^

2.*x).^2.*(-4.*sin(2.*x)./(

3.+

4.*x.^2)-32.*cos(2.*x)./(3.+4.*x.^2).^2.*x+128.*(3.^(1./2)+sin(2.*x))./(3.+4.*x.^2).^3.*x.^2.-8.*(3.^(1./2)+sin(2.*x))./(3.+4.*x.^2).^2)+16.*(1.+tan(cos((3.^(1./2)+sin(2.*x))./(3.+4.*x.^2))).^2).^2.*sin((3.^(1./2)+sin(2.*x ))./(3.+4.*x.^2)).^4.*(2.*cos(2.*x)./(3.+4.*x.^2)-8.*(3.^(1./2)+sin(2.*x))./(3.+4.*x.^2).^2.*x).^4.*tan(cos((3.^(1./2)+sin(2.*x))./(3.+4.*x.^2)))+8.*tan(cos((3.^(1./2)+sin(2.*x))./(3.+4.*x.^2))).^3.*(1.+tan(cos((3.^(1./2)+sin(2.*x))./(3.+4.*x.^2))).^2).*sin((3.^(1./2)+sin(2.*x))./(3+4.*x.^2)).^4.*(2.*cos(2.*x)./(3+4.*x.^2)-8.*(3.^(1./2)+sin(2.*x))./(3.+4.*x.^2).^2.*x).^4-8.*t an(cos((3.^(1./2)+sin(2.*x))./(3.+4.*x.^2))).*(1.+tan(cos((3.^(

1./2)+sin(

2.*x))./(

3.+4*x.^2))).^2).*sin((3.^(1./2)+sin(2.*x)).

(2.*x))./(3.+4*x.^2).^2.*x).^4+6.*(1+tan(cos((3.^(1./2)+sin(2.* x))./(3+4.*x.^2))).^2).*sin((3.^(1./2)+sin(2.*x))./(3+4.*x.^2)) .*(2.*cos(2.*x)./(3+4.*x.^2)-8.*(3.^(1./2)+sin(2.*x))./(3+4.*x. ^2).^2.*x).^2.*(-4.*sin(2.*x)./(3+4.*x.^2)-32.*cos(2.*x)./(3+4. *x.^2).^2.*x+128.*(3.^(1./2)+sin(2.*x))./(3+4.*x.^2).^3.*x.^2-8 .*(3.^(1./2)+sin(2.*x))./(3+4.*x.^2).^2)+(1+tan(cos((3.^(1./2)+ sin(2.*x))./(3+4.*x.^2))).^2).*cos((3.^(1./2)+sin(2.*x))./(3+4. *x.^2)).*(2.*cos(2.*x)./(3+4.*x.^2)-8.*(3.^(1./2)+sin(2.*x))./( 3+4.*x.^2).^2.*x).^4-3.*(1+tan(cos((3.^(1./2)+sin(2.*x))./(3+4. *x.^2))).^2).*cos((3.^(1./2)+sin(2.*x))./(3+4.*x.^2)).*(-4.*sin (2.*x)./(3+4.*x.^2)-32.*cos(2.*x)./(3+4.*x.^2).^2.*x+128.*(3.^( 1./2)+sin(2.*x))./(3+4.*x.^2).^3.*x.^2-8.*(3.^(1./2)+sin(2.*x)) ./(3+4.*x.^2).^2).^2-4.*(1+tan(cos((3.^(1./2)+sin(2.*x))./(3+4. *x.^2))).^2).*cos((3.^(1./2)+sin(2.*x))./(3+4.*x.^2)).*(2.*cos( 2.*x)./(3+4.*x.^2)-8.*(3.^(1./2)+sin(2.*x))./(3+4.*x.^2).^2.*x) .*(-8.*cos(2.*x)./(3+4.*x.^2)+96.*sin(2.*x)./(3+4.*x.^2).^2.*x+ 768.*cos(2.*x)./(3+4.*x.^2).^3.*x.^2-48.*cos(2.*x)./(3+4.*x.^2) .^2-3072.*(3.^(1./2)+sin(2.*x))./(3+4.*x.^2).^4.*x.^3+384.*(3.^

(1./2)+sin(2.*x))./(3+4.*x.^2).^3.*x)-(1+tan(cos((3.^(1./2)+sin

(2.*x))./(3+4.*x.^2))).^2).*sin((3.^(1./2)+sin(2.*x))./(3+4.*x. ^2)).*(16.*sin(2.*x)./(3+4.*x.^2)+256.*cos(2.*x)./(3+4.*x.^2).^ 2.*x-3072.*sin(2.*x)./(3+4.*x.^2).^3.*x.^2+192.*sin(2.*x)./(3+4 .*x.^2).^2-24576.*cos(2.*x)./(3+4.*x.^2).^4.*x.^3+3072.*cos(2.* x)./(3+4.*x.^2).^3.*x+98304.*(3.^(1./2)+sin(2.*x))./(3+4.*x.^2) .^5.*x.^4-18432.*(3.^(1./2)+sin(2.*x))./(3+4.*x.^2).^4.*x.^2+38 4.*(3.^(1./2)+sin(2.*x))./(3+4.*x.^2).^3)-12.*(1+tan(cos((3.^(1 ./2)+sin(2.*x))./(3+4.*x.^2))).^2).^2*sin((3.^(1./2)+sin(2.*x)) ./(3+4.*x.^2)).^2.*(2.*cos(2.*x)./(3+4.*x.^2)-8.*(3.^(1./2)+sin (2.*x))./(3+4.*x.^2).^2.*x).^4.*cos((3.^(1./2)+sin(2.*x))./(3+4 .*x.^2))-24.*tan(cos((3.^(1./2)+sin(2.*x))./(3+4.*x.^2))).^2.*( 1+tan(cos((3.^(1./2)+sin(2.*x))./(3+4.*x.^2))).^2).*sin((3.^(1. /2)+sin(2.*x))./(3+4.*x.^2)).^3.*(2.*cos(2.*x)./(3+4.*x.^2)-8.* (3.^(1./2)+sin(2.*x))./(3+4.*x.^2).^2.*x).^2.*(-4.*sin(2.*x)./( 3+4.*x.^2)-32.*cos(2.*x)./(3+4.*x.^2).^2.*x+128.*(3.^(1./2)+sin (2.*x))./(3+4.*x.^2).^3.*x.^2-8.*(3.^(1./2)+sin(2.*x))./(3+4.*x .^2).^2)-24.*tan(cos((3.^(1./2)+sin(2.*x))./(3+4.*x.^2))).^2.*( 1+tan(cos((3.^(1./2)+sin(2.*x))./(3+4.*x.^2))).^2).*sin((3.^(1. /2)+sin(2.*x))./(3+4.*x.^2)).^2.*(2.*cos(2.*x)./(3+4.*x.^2)-8.* (3.^(1./2)+sin(2.*x))./(3+4.*x.^2).^2.*x).^4.*cos((3.^(1./2)+si n(2.*x))./(3+4.*x.^2))+36.*tan(cos((3.^(1./2)+sin(2.*x))./(3+4. *x.^2))).*(1+tan(cos((3.^(1./2)+sin(2.*x))./(3+4.*x.^2))).^2).* sin((3.^(1./2)+sin(2.*x))./(3+4.*x.^2)).*(2.*cos(2.*x)./(3+4.*x .^2)-8.*(3.^(1./2)+sin(2.*x))./(3+4.*x.^2).^2.*x).^2.*cos((3.^( 1./2)+sin(2.*x))./(3+4.*x.^2)).*(-4.*sin(2.*x)./(3+4.*x.^2)-32. *cos(2.*x)./(3+4.*x.^2).^2.*x+128.*(3.^(1./2)+sin(2.*x))./(3+4. *x.^2).^3.*x.^2-8.*(3.^(1./2)+sin(2.*x))./(3+4.*x.^2).^2)+6.*ta n(cos((3.^(1./2)+sin(2.*x))./(3+4.*x.^2))).*(1+tan(cos((3.^(1./ 2)+sin(2.*x))./(3+4.*x.^2))).^2).*cos((3.^(1./2)+sin(2.*x))./(3 +4.*x.^2)).^2.*(2.*cos(2.*x)./(3+4.*x.^2)-8.*(3.^(1./2)+sin(2.* x))./(3+4.*x.^2).^2.*x).^4+6.*tan(cos((3.^(1./2)+sin(2.*x))./(3 +4.*x.^2))).*(1+tan(cos((3.^(1./2)+sin(2.*x))./(3+4.*x.^2))).^2 ).*sin((3.^(1./2)+sin(2.*x))./(3+4.*x.^2)).^2.*(-4.*sin(2.*x)./ (3+4.*x.^2)-32.*cos(2.*x)./(3+4.*x.^2).^2.*x+128.*(3.^(1./2)+si n(2.*x))./(3+4.*x.^2).^3.*x.^2-8.*(3.^(1./2)+sin(2.*x))./(3+4.* x.^2).^2).^2+8.*tan(cos((3.^(1./2)+sin(2.*x))./(3+4.*x.^2))).*( 1+tan(cos((3.^(1./2)+sin(2.*x))./(3+4.*x.^2))).^2).*sin((3.^(1. /2)+sin(2.*x))./(3+4.*x.^2)).^2.*(2.*cos(2.*x)./(3+4.*x.^2)-8.* (3.^(1./2)+sin(2.*x))./(3+4.*x.^2).^2.*x).*(-8.*cos(2.*x)./(3+4

x.^2).^3.*x.^2-48.*cos(2.*x)./(3+4.*x.^2).^2-3072.*(3.^(1./2)+s in(2.*x))./(3+4.*x.^2).^4.*x.^3+384.*(3.^(1./2)+sin(2.*x))./(3+

4.*x.^2).^3.*x)

myxx=max(yxxxx), R=(h.^4).* abs(myxx./384)

将其保存为名为myxx.m 的M 文件,然后在MATLAB 工作窗口输入该文件名

>> myxx

运行后屏幕显示myxx =)(max )4(x f x π≤≤π-和)(3,x H n 在区间],[ππ-上的误差公式

)(m a x 384

)4(4

x f h R x π≤≤π-=如下 myxx = R =

73.94706841647552 1734520780029061/9007199254740992*h^4

最后在MATLAB 工作窗口输入

>>h=2*pi/9; R =1734520780029061/9007199254740992*h^4

运行后屏幕显示)(3,x H n 在区间],[ππ-上的误差限

R =

0.04574453029948

6.7 三次样条(Spline )及其MATLAB 程序

6.7.4 用一阶导数计算的几种样条函数和MATLAB 程序

计算压紧样条的MATLAB 主程序

function

[m,H,lambda,mu,D,A,dY,sk]=ClampedSp (X,Y,dyo,dyn)

m=length(X);A=zeros(m,m);

n=m-1;H=zeros(1,n);lambda=zeros(1,n);

mu=zeros(1,n);lambda(1)=0;A(1,1)=2;A(1,2)=lambda(1);

D=zeros(1,n);H(1)=X(2)-X(1);mu(1)=1-lambda(1);D(1)=2*dyo

;

for k=1:n

hk=X(k+1)-X(k);H(k+1)=hk;

end

H=H(2:n+1);

for k=1:n-1

lambdak=H(k)/(H(k)+H(k+1)); lambda(k+1)=lambdak;

muk=1-lambda(k+1);mu(k)=muk;

dk=3*((mu(k).*(Y(k+1)-Y(k))./H(k))+(lambda(k+1).*(Y(k+2)

-Y(k+1))./H(k+1)));

D(k+1)=dk;

end

D(m)=2*dyn;mu(n)=0; n;H;lambda;mu;D;

for i=1:m-1

A(i,i)=2;A(m,m)=2; A(i,i+1)=lambda(i); A(i+1,i)=mu(i);

end

dY=A\D';

syms x

m=length(X);S=zeros(m-1,1);

for k=2:m

sk=Y(k-1)*((H(k-1)-2*X(k-1)+2*x)*(x-X(k))^2)/(H(k-1)^3)+

Y(k)*((H(k-1)+2*X(k)-2*x)*(x-X(k-1))^2)/(H(k-1)^3)+dY(k-1)*((x-X(k-1))*(x-X(k))^2)/(H(k-1)^2)+dY(k)*((x-X(k))*(x-X(k-1))^2)/(H (k-1)^2)

end

例6.7.2 用计算压紧样条的MATLAB 主程序ClampedSp.m 求例6.7.1的问题.

解 保存ClampedSp.m 为M 文件,输入程序

>> X=[0:2:6,10],Y=[0,16,36,54,82],dyo=8,dyn=7,

[m,H,lambda,mu,D,A,dY,sk]=ClampedSp(X,Y,dyo,dyn)

运行后输出压紧样条函数sk ,X 的维数m ,由1-i h ),,2,1(n i =、i λ、i μ和i d

)1,,2,1(-=n i 的坐标组成的向量H 、lambda 和mu 、D ,

(6.90)式的系数矩阵A ,线性方程组(6.90)的解向量dY 如下

sk =

2*(6-2*x)*x^2+2*x*(x-2)^2+9/4*(x-2)*x^2

sk =

2*(-2+2*x)*(x-4)^2+9/2*(10-2*x)*(x-2)^2+9/4*(x-2)*(x-4)^

2+5/2*(x-4)*(x-2)^2

sk =

9/2*(-6+2*x)*(x-6)^2+27/4*(14-2*x)*(x-4)^2+5/2*(x-4)*(x-

6)^2+2*(x-6)*(x-4)^2

sk =

27/32*(-8+2*x)*(x-10)^2+41/32*(24-2*x)*(x-6)^2+1/2*(x-6)

*(x-10)^2+7/16*(x-10)*(x-6)^2

m =

5

H =

2 2 2 4

lambda =

0 0.5000 0.5000 0.3333

mu =

0.5000 0.5000 0.6667 0

D =

16.0000 27.0000 28.5000 25.0000 14.0000

A =

2.0000 0 0 0 0

0.5000 2.0000 0.5000 0 0

0 0.5000 2.0000 0.5000 0

0 0 0.6667 2.0000 0.3333

0 0 0 0 2.0000

dY =

8.0000

9.0000

10.0000

8.0000

7.0000

输入程序

>> x2=3,

s2=2*(-2+2*x)*(x-4)^2+9/2*(10-2*x)*(x-2)^2+9/4*(x-2)*(x-

4)^2+5/2*(x-4)*(x-2)^2

x4=8,

s4=27/32*(-8+2*x)*(x-10)^2+41/32*(24-2*x)*(x-6)^2+1/2*(x

-6)*(x-10)^2+7/16*(x-10)*(x-6)^2

运行后输出)3()3(2s f ≈和)8()8(4s f ≈的值如下

x2 = s2 = x4 = s4 =

3 25.7500 8 68.5000

例6.7.2和例6.7.1计算的结果完全相同.

计算自然样条的MATLAB 主程序

function [m,H,lambda,mu,D,A,dY,sk]=NaturalSp(X,Y)

m=length(X);A=zeros(m,m);

n=m-1;H=zeros(1,n);lambda=zeros(1,n);

mu=zeros(1,n);lambda(1)=1;A(1,1)=2;A(1,2)=lambda(1);

D=zeros(1,n);H(1)=X(2)-X(1);mu(1)=1;D(1)=3*(Y(2)-Y(1));

for k=1:n

end

H=H(2:n+1);

for k=1:n-1

lambdak=H(k)/(H(k)+H(k+1)); lambda(k+1)=lambdak;

muk=1-lambda(k+1);mu(k)=muk;

dk=3*((mu(k).*(Y(k+1)-Y(k))./H(k))+(lambda(k+1).*(Y(k+2)

-Y(k+1))./H(k+1)));

D(k+1)=dk;

end

D(m)= 3*(Y(m)-Y(m-1))/H(m-1);mu(n)=1; n;H;lambda;mu;D;

for i=1:m-1

A(i,i)=2;A(m,m)=2; A(i,i+1)=lambda(i); A(i+1,i)=mu(i);

end

dY=A\D';

syms x

m=length(X);S=zeros(m-1,1);

for k=2:m

sk=Y(k-1)*((H(k-1)-2*X(k-1)+2*x)*(x-X(k))^2)/(H(k-1)^3)+

Y(k)*((H(k-1)+2*X(k)-2*x)*(x-X(k-1))^2)/(H(k-1)^3)+dY(k-1)*((x-X(k-1))*(x-X(k))^2)/(H(k-1)^2)+dY(k)*((x-X(k))*(x-X(k-1))^2)/(H (k-1)^2)

end

例6.7.3 用计算自然样条的MATLAB 主程序NaturalSp.m 求例6.7.1的问题.

解 保存ClampedSp.m 为M 文件,输入程序

>> X=[0:2:6,10],Y=[0,16,36,54,82],dyo=8,dyn=7,

[m,H,lambda,mu,D,A,dY,sk]=NaturalSp(X,Y)

运行后输出自然样条函数sk , X 的维数m ,由1-i h ),,2,1(n i =、i λ、i μ和i d )1,,2,1(-=n i 的坐标组成的向量H 、lambda 和mu 、D ,(6.88)式的系数矩阵A ,线性方程组(6.88)的解向量dY 如下

sk =

2*(6-2*x)*x^2+915/172*x*(x-2)^2+117/86*(x-2)*x^2

sk =

2*(-2+2*x)*(x-4)^2+9/2*(10-2*x)*(x-2)^2+117/86*(x-2)*(x-

4)^2+471/172*(x-4)*(x-2)^2

sk =

9/2*(-6+2*x)*(x-6)^2+27/4*(14-2*x)*(x-4)^2+471/172*(x-4)

*(x-6)^2+333/172*(x-6)*(x-4)^2

sk =

27/32*(-8+2*x)*(x-10)^2+41/32*(24-2*x)*(x-6)^2+333/688*(

x-6)*(x-10)^2+285/688*(x-10)*(x-6)^2

m = 5

H = 2 2 2 4

lambda = 1 1/2 1/2 1/3

mu = 1/2 1/2 2/3 1

D = 48 27 57/2 25 21

A =

2 1 0 0 0

1/2 2 1/2 0 0

0 1/2 2 1/2 0

0 0 2/3 2 1/3

0 0 0 1 2

dY =

915/43

234/43

471/43

333/43

计算方法 课内实验 插值法与函数逼近

《计算方法》课内实验报告 学生姓名:张学阳1009300132 及学号: 学院: 理学院 班级: 数学101 课程名称:计算方法 实验题目:插值法与函数逼近 指导教师 宋云飞讲师 姓名及职称: 朱秀丽讲师 尚宝欣讲师 2012年10月15日

目录 一、实验题目.......................................................... 错误!未定义书签。 二、实验目的.......................................................... 错误!未定义书签。 三、实验内容.......................................................... 错误!未定义书签。 四、实现结果.......................................................... 错误!未定义书签。 五、实验体会或遇到问题 (6)

插值法与函数逼近 二、实验目的 1.熟悉matlab 编写及运行数值计算程序的方法。 2.进一步理解插值法及函数逼近方法的理论基础。 3.进一步掌握给定数据后应用插值法及函数逼近方法进行数据处理并给出图示结果的实际操作过程。 三、实验内容 1.已知函数在下列各点的值为 试用4次牛顿插值多项式)(4x P 及三次样条函数)(x S (自然边界条件)对数据进行插值。给出求解过程,并用图给出 (){},10,1,0),()(,08.02.0,,4 ===+=i x S y x P y i x y x i i i i i 及。 2.下列数据点的插值 可以得到平方根函数的近似。 (1)用这9个点作8次多项式插值)(8x L 。 (2)用三次样条(第一类边界条件)插值给出)(x S 。 给出求解过程,在区间[0,64]上作图,从得到的结果看,在区间[0,64]上哪种插值结果更精确?在区间[0,1]上两种插值哪个更精确? 3.由实验给出数据表 试求3次、4次多项式的曲线拟合,再根据数据曲线形状,求一个另外函数的拟合曲线。给出求解过程,用图表示实验数据曲线及三种拟合曲线。

求多元函数极限的方法

求多元函数极限的方法 【摘要】对于大部分学生,尤其是初接触高等数学的同学而言,极限是一道很难过的关,因为那种“无限逼近”却又“无法达到”的抽象对于刚刚结束中学数学学习,习惯于具体图形分析、函数计算的同学来说,在思维上有了更高的要求。而对于高等数学来讲,极限又是相当重要的基础,不管是函数连续性的验证,亦或是单侧导数的求解,极限都是很重要的一个环节,它就相当于一条线惯于始终,所以说学好极限,是学好高等数学的一个起点。【1】 【关键词】多元函数;求极限多种方法;求极限常出现的错误 【引言】之前学过如连续、导数微分和积分等都要用极和秋极限的方法,例如:利用定义来求极限、用柯西收敛准则、利用两边夹定理等等。这些方法虽然简便易于理解和掌握,但对 于一些特殊的极限题目很难解决,例如:设0a >,10a >,2 12(3) 3n n n n a a a a a a ++=+求lim n n a →∞的问题题目尽给出了第n 项和第n +1项的关系若用利用定义来求极限、用柯西收敛准则 1 ! lim ! n k n k n =→∞ ∑及求一些复合函数极限的问题本文将探讨一些特殊的求极限的方法,对某些用常 见方法不易求解的题目运用此方法可以容易地解出。【2】本文将从多个方面,通过利用极限的性质及相关概念和几个典型例题对常用求极限的方法进行解析,并列出容易出错的地方。 1 利用极限定义的思想观察函数的极限 例1、讨论当x → 12时函数y =21 x x +的极限。我们列出了当x →12 时某些函数值,考察 从列表可以看出,当x 趋向于2时,y 就趋向于0.7,即x →2 时,y =21 x x +的极限是0.75。 2、利用四则运算法则求极限 例2(1)求2 3 32 1 lim(4)x x x →-+ (2)221 lim 21 x x x →-+ 解(2)2 21lim 21x x x →-+=2 2 2 lim(1)3lim(21)5 x x x x →→-=+ 3、利用无穷小量与无穷大量的关系及无穷小量的性质求极限 例3求0 1 lim sin x x x →

实验四插值法

实验四、插值法 插值法是函数逼近的一种重要方法,它是数值积分、微分方程数值解等数值计算的基础与工具,其中多项式插值是最常用和最基本的方法。拉格朗日插值多项式的优点是表达式简单明确,形式对称,便于记忆,它的缺点是如果想要增加插值节点,公式必须整个改变,这就增加了计算工作量。而牛顿插值多项式对此做了改进,当增加一个节点时只需在原牛顿插值多项式基础上增加一项,此时原有的项无需改变,从而达到节省计算次数、节约存储单元、应用较少节点达到应有精度的目的。 一、实验目的 1、理解插值的基本概念,掌握各种插值方法,包括拉格朗日插值和牛顿插值等,注意其不同特点; 2、通过实验进一步理解并掌握各种插值的基本算法。 二、Matlab命令和程序 命令poly:创建一个向量,其分量为一个多项式的系数,该多项式具有给定的根。 命令polyval:求多项式的值, 命令 conv: 创建一个向量,其分量为一个多项式的系数,该多项式是另外两个多项式的积 polyval(C,2> >> P=poly(2> P=1 -2

Q=poly(3> Q=1 -3 >> conv(P,Q> ans= 1 -5 6 >> polyval(P,2> ans= 1、拉格朗日插值( 基于N+1个点,计算拉格朗日多项式> function [C,L]=lagran(X,Y> %input --X is a vector that contains a list of abscissasb5E2RGbCAP % Y is a vector that contains a list of ordinatesp1EanqFDPw %output--C is a matrix that contains the coefficient of the lagraneDXDiTa9E3d % interplatory polynomial % -- L is a matrix that contains the Lagrange coefficent polynomialsRTCrpUDGiT w=length(X>。 n=w-1。

MATLAB数值实验一(数据的插值运算及其应用完整版)

佛山科学技术学院 实 验 报 告 课程名称 数值分析 实验项目 插值法与数据拟合 专业班级 机械工程 姓 名 余红杰 学 号 10 指导教师 陈剑 成 绩 日 期 月 日 一、实验目的 1、学会Lagrange 插值、牛顿插值和三次样条插值等基本插值方法; 2、讨论插值的Runge 现象 3、学会Matlab 提供的插值函数的使用方法,会用这些函数解决实际问题。 二、实验原理 1、拉格朗日插值多项式 2、牛顿插值多项式 3、三次样条插值 三、实验步骤 1、用MATLAB 编写独立的拉格朗日插值多项式函数 2、用MATLAB 编写独立的牛顿插值多项式函数 3、用MATLAB 编写独立的三次样条函数(边界条件为第一、二种情形) 4、已知函数在下列各点的值为: 根据步骤1,2,3编好的程序,试分别用4次拉格朗日多项式4()L x 、牛顿插值多项式4()P x 以及三次样条函数()S x (自然边界条件)对数据进行插值,并用图给出 {(,),0.20.08,0,1,2, ,10i i i x y x i i =+=},4()L x 、4()P x 和()S x 。 5、在区间[-1,1]上分别取10,20n =用两组等距节点对龙格函数 2 1 (),(11)125f x x x = -≤≤+作多项式插值,对不同n 值,分别画出插值函数及()f x 的图形。 6、下列数据点的插值

可以得到平方根函数的近似,在区间[0,64]上作图。 (1)用这9个点作8次多项式插值8()L x 。 (2)用三次样条(第一边界条件)程序求()S x 。 7、对于给函数2 1 ()125f x x = +在区间[-1,1]上取10.2(0,1, ,10)i x i i =-+=,试求3次 曲线拟合,试画出拟合曲线并打印出方程,与第5题的结果比较。 四、实验过程与结果: 1、Lagrange 插值多项式源代码: function ya=lag(x,y,xa) %x 所有已知插值点 %y 插值点对应函数值 %xa 所求点,自变量 %ya 所求点插值估计量 ya=0; mu=1; %初始化 %循环方式求L 系数,并求和: for i = 1:length(y) for j = 1:length(x) if i ~= j mu = mu * (xa - x(j) ) / ( x(i) - x(j) ); else continue end end ya = ya + y(i) * mu ; mu = 1; end 2、Newton 源代码: function ya = newton(x,y,xa) %x 所有已知插值点 %y 插值点对应函数值 %xa 所求点,自变量 %ya 所求点插值估计量 %建立系数零矩阵D 及初始化:

插值法实验报告

实验二插值法 1、实验目的: 1、掌握直接利用拉格郎日插值多项式计算函数在已知点的函数值;观察拉格郎日插值的龙格现象。 2、了解Hermite插值法、三次样条插值法原理,结合计算公式,确定函数值。 2、实验要求: 1)认真分析题目的条件和要求,复习相关的理论知识,选择适当的解决方案和算法; 2)编写上机实验程序,作好上机前的准备工作; 3)上机调试程序,并试算各种方案,记录计算的结果(包括必要的中间结果); 4)分析和解释计算结果; 5)按照要求书写实验报告; 3、实验内容: 1) 用拉格郎日插值公式确定函数值;对函数f(x)进行拉格郎日插值,并对f(x)与插值多项式的曲线作比较。 已知函数表:(0.56160,0.82741)、(0.56280,0.82659)、(0.56401,0.82577)、(0.56521,0.82495)用三次拉格朗日插值多项式求x=0.5635时函数近似值。 2) 求满足插值条件的插值多项式及余项 1) 4、题目:插值法 5、原理: 拉格郎日插值原理: n次拉格朗日插值多项式为:L n (x)=y l (x)+y 1 l 1 (x)+y 2 l 2 (x)+…+y n l n (x)

n=1时,称为线性插值, L 1(x)=y (x-x 1 )/(x -x 1 )+y 1 (x-x )/(x 1 -x )=y +(y 1 -x )(x-x )/(x 1 -x ) n=2时,称为二次插值或抛物线插值, L 2(x)=y (x-x 1 )(x-x 2 )/(x -x 1 )/(x -x 2 )+y 1 (x-x )(x-x 2 )/(x 1 -x )/(x 1 -x 2 )+y 2 (x -x 0)(x-x 1 )/(x 2 -x )/(x 2 -x 1 ) n=i时, Li= (X-X0)……(X-X i-1)(x-x i+1) ……(x-x n) (X-X0)……(X-X i-1)(x-x i+1) ……(x-x n) 6、设计思想: 拉格朗日插值法是根据n + 1个点x0, x1, ... x n(x0 < x1 < ... x n)的函数值f (x0), f (x1) , ... , f (x n)推出n次多項式p(x),然后n次多項式p (x)求出任意的点x对应的函数值f (x)的算法。 7、对应程序: 1 ) 三次拉格朗日插值多项式求x=0.5635时函数近似值 #include"stdio.h" #define n 5 void main() { int i,j; float x[n],y[n]; float x1; float a=1; float b=1; float lx=0; printf("\n请输入想要求解的X:\n x="); scanf("%f",&x1); printf("请输入所有点的横纵坐标:\n"); for(i=1;i

Matlab中插值函数汇总和使用说明.

告: Matlab中插值函数汇总和使用说明收藏 命令1 interp1 功能一维数据插值(表格查找。该命令对数据点之间计算内插值。它找出一元函数f(x在中间点的数值。其中函数f(x由所给数据决定。x:原始数据点 Y:原始数据点 xi:插值点 Yi:插值点 格式 (1yi = interp1(x,Y,xi 返回插值向量yi,每一元素对应于参量xi,同时由向量x 与Y 的内插值决定。参量x 指定数据Y 的点。 若Y 为一矩阵,则按Y 的每列计算。yi 是阶数为length(xi*size(Y,2的输出矩阵。 (2yi = interp1(Y,xi 假定x=1:N,其中N 为向量Y 的长度,或者为矩阵Y 的行数。 (3yi = interp1(x,Y,xi,method 用指定的算法计算插值: ’nearest’:最近邻点插值,直接完成计算; ’linear’:线性插值(缺省方式,直接完成计算;

’spline’:三次样条函数插值。对于该方法,命令interp1 调用函数spline、ppval、mkpp、umkpp。这些命令生成一系列用于分段多项式操作的函 数。命令spline 用它们执行三次样条函数插值; ’pchip’:分段三次Hermite 插值。对于该方法,命令interp1 调用函数p chip,用于对向量x 与y 执行分段三次内插值。该方法保留单调性与数据的外形; ’cubic’:与’pchip’操作相同; ’v5cubic’:在MATLAB 5.0 中的三次插值。 对于超出x 范围的xi 的分量,使用方法’nearest’、’linear’、’v5cubic’的插值算法,相应地将返回NaN。对其他的方法,interp1 将对超出的分量执行外插值算法。 (4yi = interp1(x,Y,xi,method,'extrap' 对于超出x 范围的xi 中的分量将执行特殊的外插值法extrap。 (5yi = interp1(x,Y,xi,method,extrapval 确定超出x 范围的xi 中的分量的外插值extrapval,其值通常取NaN 或0。 例1 1.>>x = 0:10; y = x.*sin(x; 2.>>xx = 0:.25:10; yy = interp1(x,y,xx; 3.>>plot(x,y,'kd',xx,yy 复制代码 例2 1.>> year = 1900:10:2010;

求二元函数极限地几种方法

精彩文档 1.二元函数极限概念分析 定义1 设函数f 在2D R ?上有定义,0P 是D 的聚点,A 是一个确定的实数.如果对于任意给定的正数ε,总存在某正数δ,使得00(;)P U P D δ∈时,都有 ()f P A ε-<, 则称f 在D 上当0P P →时,以A 为极限,记0 lim ()P P P D f P A →∈=. 上述极限又称为二重极限. 2.二元函数极限的求法 2.1 利用二元函数的连续性 命题 若函数(,)f x y 在点00(,)x y 处连续,则 0000(,)(,) lim (,)(,)x y x y f x y f x y →=. 例1 求2 (,)2f x y x xy =+ 在点(1,2)的极限. 解: 因为2 (,)2f x y x xy =+在点(1,2)处连续,所以 12 212 2lim (,) lim(2) 12125.x y x y f x y x xy →→→→=+=+??= 例2 求极限()()2 21,1,21 lim y x y x +→. 解: 因函数在()1,1点的邻域内连续,故可直接代入求极限,即 ()()221,1,21 lim y x y x +→=3 1.

精彩文档 2.2 利用恒等变形法 将二元函数进行恒等变形,例如分母或分子有理化等. 例3 求 00 x y →→ 解: 00 x y →→ 00 x y →→= 00 x y →→= 00 1. 4 x y →→==-例4 ()() 2 2 220,0,321 )31)(21(lim y x y x y x +-++→. 解: 原式()() ( ) ) () () ,0,02 211lim 231x y x y →+= ++ ()( 22 ,0,0lim x y →= + 11022 = +=.

计算方法实验

算方法实验指导 姓名学号院系专业哈尔滨工业大学

计算方法实验指导 根据实际问题建立的数学模型,一般不能求出所谓的解析解,必须针对数学模型 的特点确定适当的计算方法,编制出计算机能够执行的计算程序,输入计算机,进行 调试,完成运算,如果计算结果存在问题或不知是否正确,还需要重新确定新的计算 方法,再编制出计算程序,输入计算机,重新调试,完成运算,直至获得正确的计算 结果,这就是数值计算的全部过程。 学生在学习“计算方法”和“高级语言”等课程时普遍存在的问题是:只会套用 教科书中的标准程序进行数值计算,很少有人能够独立地将学过的数值算法编制成计 算机程序,至于灵活应用已经掌握的算法求解综合性较大的课题,则更是困难的事情。 编写《计算方法实验指导》的目的是:突出数值计算程序结构化的思想。提高学 生的编程能力,加深对“计算方法”课程内容的理解和掌握,为”计算方法“课程的 教学服务,进一步奠定从事数值计算工作的基础。具体地 1. 根据“计算方法”课程内容的特点,给出五个典型算法的分析流程,学生可以 利用所掌握的 “高级语言”顺利地编制出计算机程序,上机实习,完成实验环节的教 学要求。 2. 所有的计算实习题目都经过任课教师逐一检验,准确无误。 3. 充分利用循环的思想、 迭代的思想, 给出算法结构描述和程序语言的对应关系, 有利于学生编 制相应的程序。 4. 结合实习题目,提出实验要求,要求学生按规范格式写出相应的实验报告,实 验报告成绩记入 期末总成绩。需要提醒学生:不能简单地套用现成的标准程序完成实 验题目,应当把重点放在对算法的理解、程序的优化设计、上机调试和计算结果分析 上,否则就失去实验课的目的啦。 5. 五个具体的实验题目是: 实验题目 实验题目 实验题目 实验题目 实验题目 要求必须完 成其中三个(如果全部完成更好) 。 1 拉格朗日 (Lagrange) 插值 2 龙贝格 (Romberg) 积分法 3 四阶龙格—库塔 (Runge — Kutta) 方法 4 牛顿 (Newton) 迭代法 5 高斯 (Gauss) 列主元消去法

数值分析(计算方法)实验一

《数值分析》 课程实验指导书 实验一 函数插值方法 一、问题提出 对于给定的一元函数)(x f y =的n+1个节点值(),0,1,,j j y f x j n == 。试用Lagrange 公式求其插值多项式或分段二次Lagrange 插值多项式。 数据如下: (1) j x 0.4 0.55 0.65 0.80 0.95 1.05 j y 0.41075 0.57815 0.69675 0.90 1.00 1.25382 求五次Lagrange 多项式5L ()x ,和分段三次插值多项式,计算(0.596)f ,(0.99)f 的值。(提示:结果为(0.596)0.625732f ≈, (0.99) 1.05423f ≈ ) (2) j x 1 2 3 4 5 6 7 j y 0.368 0.135 0.050 0.018 0.007 0.002 0.001 试构造Lagrange 多项式6L ()x ,计算的(1.8)f ,(6.15)f 值。(提示:结果为(1.8)0.164762f ≈, (6.15)0.001266f ≈ ) 二、要求 1、 利用Lagrange 插值公式 00,()n n i n k k i i k k i x x L x y x x ==≠??-= ?-??∑∏编写出插值多项式程序; 2、 给出插值多项式或分段三次插值多项式的表达式; 3、 根据节点选取原则,对问题(2)用三点插值或二点插值,其结果如何; 4、 对此插值问题用Newton 插值多项式其结果如何。

四、实验分析: Lagrange 插值多项式的表达式: 1,,2,1,)()()(, )()(1111+=--==∏∑+≠=+=n i x x x x x l x l y x L n i j j j i j i n i i i 。 其中)(x l i 被称为插值基函数,实际上是一个n 次多项式。)(x l i 的这种表示具有较好的对称性。公式具有两大优点:(1)求插值多项式,不需要求解线性方程组,当已知数据点较多时,此公式更能显示出优越性。(2)函数值可以用符号形式表示,数据点未确定的纵坐标可用多项式表示。 Newton 插值多项式如下: 10010,()()[,,]()k n n j k k j j k N x f x f x x x x -==≠=+?-∑∏ 其中: 00,0()()[,,]k i k i i j j j i k f x x x f x x ==≠-=∑∏ Newton 插值多项式的优点是:当每增加一个节点时,只增加一项多项式。 三、实验程序及注释 1、m 程序: function [c,l]=lagran(x,y) % x 为n 个节点的横坐标组成的向量,y 为纵坐标所组成的向量 % c 为所得插值函数的系数所组成的向量 w=length(x); n=w-1; l=zeros(w,w); for k=1:n+1 v=1; for j=1:n+1 if k~=j v=conv(v,poly(x(j)))/(x(k)-x(j)); end end l(k,:)=v; end c=y*l; function fi=Lagran_(x,f,xi) fi=zeros(size(xi)); n=length(f); for i=1:n

[整理]matlab绘制二元函数图形.

MATL AB绘制二元函数的图形 【实验目的】 1.了解二元函数图形的绘制。 2.了解空间曲面等高线的绘制。 3.了解多元函数插值的方法。 4.学习、掌握MATLAB软件有关的命令。 【实验内容】 画出函数2 2y z+ =的图形,并画出其等高线。 x 【实验准备】 1.曲线绘图的MATLAB命令 MATLAB中主要用mesh,surf命令绘制二元函数图形。主要命令mesh(x,y,z)画网格曲面,这里x,y,z是数据矩阵,分别表示数据点的横坐标,纵坐标和函数值,该命令将数据点在空间中描出,并连成网格。 surf(x,y,z)画完整曲面,这里x,y,z是数据矩阵,分别表示数据点的横坐标,纵坐标和函数值,该命令将数据点所表示曲面画出。 【实验重点】 1. 二元函数图形的描点法 2. 曲面交线的计算 3. 地形图的生成

【实验难点】 1. 二元函数图形的描点法 2. 曲面交线的计算 【实验方法与步骤】 练习1画出函数2 2y =的图形,其中]3,3 x z+ ? - y x。 ∈ , [ ]3,3 [ (- ) 用MATLAB作图的程序代码为 >>clear; >>x=-3:0.1:3; %x的范围为[-3,3] >>y=-3:0.1:3; %y的范围为[-3,3] >>[X,Y]=meshgrid(x,y); %将向量x,y指定的区域转化为矩阵X,Y >>Z=sqrt(X.^2+Y.^2); %产生函数值Z >>mesh(X,Y,Z) 运行结果为

图5.3 如果画等高线,用contour,contour3命令。 contour画二维等高线。 contour3画三维等高线。画图5.3所示的三维等高线的MA TLAB 代码为 >>clear; >>x=-3:0.1:3; >>y=-3:0.1:3; >>[X,Y]=meshgrid(x,y); >>Z=sqrt(X.^2+Y.^2); >>contour3(X,Y,Z,10); %画10条等高线 >>xlabel('X-axis'),ylabel('Y-axis'),zlabel('Z-axis'); %三个坐标轴的

数学分析下——二元函数的极限课后习题

第二节 二元函数的极限 1、试求下列极限(包括非正常极限): (1)(,)(0,0)lim x y x 2y 2x 2+y 2 ; (2)(,)(0,0)lim x y 1+x 2+y 2x 2+y 2 ; (3) (,)(0,0) lim x y x 2+y 21+x 2+y 2 -1 ; (4)(,)(0,0)lim x y xy+1 x 4+y 4 ; (5)(,) (1,2)lim x y 12x-y ; (6)(,)(0,0)lim x y (x+y)sin 1 x 2+y 2 ; (7)(,)(0,0)lim x y sin(x 2+y 2)x 2+y 2 x 2+y 2 . 2、讨论下列函数在点(0,0)的重极限与累次极限: (1)f(x,y)=y 2x 2+y 2 ; (2)f(x,y)=(x+y)sin 1x sin 1y ; (3)f(x,y)=x 2y 2x 2y 2+(x-y)2 ; (4)f(x,y)=x 3+y 3 x 2+y ; (5)f(x,y)=ysin 1x ; (6)f(x,y)=x 2y 2 x 3+y 3 ; (7)f(x,y)=e x -e y sinxy . 3、证明:若1 。 (a,b) lim (x,y )f(x,y)存在且等于A ;2。 y 在b 的某邻域内,有 lim x a f(x,y)= (y)则 y b lim a lim x f(x,y)=A. 4、试应用ε—δ定义证明 (x,y)(0,0)lim x 2y x 2+y 2 =0. 5、叙述并证明:二元函数极限的唯一性定理、局部有界性定理与局部保号性定理. 6、试写出下列类型极限的精确定义: (1) (x,y) ( ,) lim f(x,y)=A ; (2) (x,y) (0, ) lim f(x,y)=A. 7、试求下列极限: (1) (x,y) ( , )lim x 2+y 2 x 4+y 4 ; (2)(x,y)(, ) lim (x 2+y 2)e -(x+y);

计算方法-插值方法实验

实验一插值方法 一. 实验目的 (1)熟悉数值插值方法的基本思想,解决某些实际插值问题,加深对数值插值方法 的理解。 (2)熟悉Matlab 编程环境,利用Matlab 实现具体的插值算法,并进行可视化显示。 二. 实验要求 用Matlab 软件实现Lagrange 插值、分段线性插值、三次Hermite 插值、Aitken 逐步插值算法,并用实例在计算机上计算和作图。 三. 实验内容 1. 实验题目 (1 ) 已 知概 率积 分dx e y x x ?-= 2 2 π 的数据表 构造适合该数据表的一次、二次和三次Lagrange 插值公式,输出公式及其图形,并计算x =0.472时的积分值。 答: ①一次插值公式: 输入下面内容就可以得到一次插值结果 >> X=[0.47,0.48];Y=[0.4937452,0.5027498]; >> x=0.472; >> (x-X(2))/(X(1)-X(2))*Y(1)+(x-X(1))/(X(2)-X(1))*Y(2) ans =0.495546120000000 >> ②两次插值公式为: 输入下面内容就可以得到两次插值结果 >> X=[0.46,0.47,0.48];Y=[0.4846555,0.4937452,0.5027498]; >> x=0.472; >>(x-X(2))*(x-X(3))/((X(1)-X(2))*(X(1)-X(3)))*Y(1)+(x-X(1))*(x-X(3))/((X(2)-X(1))*(X(2)-X(3)))*Y(2)+(x-X(2))*(x-X(1))/((X(3)-X(2))*(X(3)-X(1)))*Y(3) i 0 1 2 3 x 0.46 047 0.48 0.49 y 0.4846555 0.4937452 0.5027498 0.5116683

(整理)二元函数极限的求法.

二元函数极限的求法 数学与统计学院、数学与应用数学、0701班,湖北,黄石,435002 1.引言 多元函数的极限在高等数学中非常重要,但由于多元函数的自变量多,因此对于判断其极限存在与否及其求法,比起一元函数的极限就显得比较困难.求极限和证明极限的方法很多,一般我们常用定义法,初等变形法,两边夹准则,阶的估计等.在这几种方法中,定义法是基础,但是比较繁琐,其他方法有的较易,有的较难,让人不知道从何下手.因此,我们有必要总结探讨出比较容易好的方法去求多元函数的极限.多元函数极限在现在的生活中也有很大的用处,比如工程计算方面.从以上来看,研究归纳总结多元函数极限的求法问题是有意义和必要的.本文主要研究二元函数极限的定义以及二元函数极限求解的几种方法,并以实例加以说明. 2.二元函数极限的定义 定义1 设E 是2R 的一个子集,R 是实数集,f 是一个规律,如果对E 中的每一点(,)x y ,通过规律f ,在R 中有唯一的一个u 与此对应,则称f 是定义在E 上的一个二元函数,它在点(,)x y 的函数值是u ,并记此值为(,)f x y ,即(,)u f x y =. 有时,二元函数可以用空间的一块曲面表示出来,这为研究问题提供了直观想象.例如,二元函数222y x R x --=就是一个上半球面,球心在原点,半径为R ,此函数定义域为满足关系式222R y x ≤+的x ,y 全体,即 }|),{(222R y x y x D ≤+=.又如,xy Z =是马鞍面. 知道多元函数的定义之后,在我们求多元函数极限之前我们必须知道多

元函数极限的定义. 定义2 设E 是2R 的一个开集,A 是一个常数,二元函数()(,)f M f x y =在点()000,M x y E ∈附近有定义.如果0>?ε,0>?δ,当()00,r M M δ<<时,有()f M A ε-<,就称A 是二元函数在0M 点的极限.记为()0 lim M M f M A →=或 ()()0f M A M M →→. 定义的等价叙述 1 :设E 是2R 的一个开集,A 是一个常数,二元函数 ()(,)f M f x y =在点()000,M x y E ∈附近有定义.如果0>?ε,0>?δ,当()() 22 000x x y y δ< -+-<时,有(,)f x y A ε-<,就称A 是二元函数在0 M 点的极限。记为()0 lim M M f M A →=或()()0f M A M M →→. 定义的等价叙述2: 设E 是2R 的一个开集,A 是一个常数,二元函数 ()(,)f M f x y =在点()000,M x y E ∈附近有定义.如果0>?ε,0>?δ,当 000,0x x y y δδ<-<<-<且()()00,,x y x y ≠时, 有(,)f x y A ε-<,就称A 是二元函数在0M 点的极限.记为 ()0 l i m M M f M A →=或 ()()0f M A M M →→. 注:(1)和一元函数的情形一样,如果0 lim ()M M f M A →=,则当M 以任何 点列及任何方式趋于0M 时,()f M 的极限是A ;反之,M 以任何方式及任何点列趋于0M 时,()f M 的极限是A .但若M 在某一点列或沿某一曲线0M →时,()f M 的极限为A ,还不能肯定()f M 在0M 的极限是A . 二元函数的极限较之一元函数的极限而言,要复杂得多,特别是自变量的变化趋势,较之一元函数要复杂.

插值与多项式逼近的数组计算方法实验讲解

插值与多项式逼近的数组计算方法实验 郑发进 2012042020022 【摘要】计算机软件中经常要用到库函数,如) cos,x e,它们 (x (x sin,) 是用多项式逼近来计算的。虽然目前最先进的逼近方法是有理函数(即多项式的商),但多项式逼近理论更适于作为数值分析的入门课程。在已知数据具有高精度的情况下,通常用组合多项式来构造过给定数据点的多项式。构造组合多项式的方法有许多种,如线性方程求解、拉格朗日系数多项式以及构造牛顿多项式的方分和系数表。 关键字泰勒级数、拉格朗日插值法、牛顿插值法、帕德逼近 一、实验目的 1.通过具体实验,掌握泰勒级数、拉格朗日插值法、牛顿插值法、帕德逼近的编程技巧。 2.比较各插值方法的优劣并掌握。 二、实验原理 1.泰勒级数 在数学中,泰勒级数(英语:Taylor series)用无限项连加式——级数来表示一个函数,这些相加的项由函数在某一点的导数求得。 如果在点x=x 具有任意阶导数,则幂级数 称为在点x 处的泰勒级数。 =0,得到的级数 在泰勒公式中,取x 称为麦克劳林级数。函数的麦克劳林级数是x的幂级数,那么这种展开

是唯一的,且必然与的麦克劳林级数一致。 2.拉格朗日插值法 如对实践中的某个物理量进行观测,在若干个不同的地方得到相应的观测值,拉格朗日插值法可以找到一个多项式,其恰好在各个观测的点取到观测到的值。这样的多项式称为拉格朗日(插值)多项式。数学上来说,拉格朗日插值法可以给出一个恰好穿过二维平面上若干个已知点的多项式函数。 在平面上有(x 1,y 1)(x 2,y 2)...(x n ,y n )共n 个点,现作一条函数f (x )使其图像经过这n 个点。 作n 个多项式p i (x),i=1,2,3...,n,使得 最后可得 3.牛顿插值法 插值法利用函数f (x)在某区间中若干点的函数值,作出适当的特定函数,在这些点上取已知值,在区间的其他点上用这特定函数的值作为函数f (x)的近似值。如果这特定函数是多项式,就称它为插值多项式。利用插值基函数很容易得到拉格朗日插值多项式,公式结构紧凑,在理论分析中甚为方便,但当插值节点增减时全部插值基函数均要随之变化,整个公式也将发生变化, 这在实际计算中是很不方便的,为了克服这一缺点,提出了牛顿插值。 牛顿插值通过求各阶差商,递推得到的一个公式: 10121()()()()()()N N N N P x P x a x x x x x x x x --=+---- 牛顿插值与拉格朗日插值具有唯一性。 4.帕德逼近 它不仅与逼近论中其他许多方法有着密切的关系,而且在实际问题特别是许多物理问题中有着广泛的应用。设是在原点某邻域内收敛的、具有复系数的麦克劳林级数。欲确定一个有理函数,式中,使得前次方的系数为0,即使得 此处约定qk =0(k>n )。虽然所求得的Pm(z)和Qn(z)不惟一,但是比式却总是惟一的。有理函数称为F(z)的(m,n)级帕德逼近,记为(m/n)。由(m/n)所形成的阵列称为帕德表。

实验5 插值方法

实验5 插值方法 一、实验目的及意义 [1] 了解插值的基本原理 [2] 了解拉格朗日插值、线性插值、样条插值的基本思想; [3] 了解三种网格节点数据的插值方法的基本思想; [4] 掌握用MATLAB 计算三种一维插值和两种二维插值的方法; [5] 通过范例展现求解实际问题的初步建模过程; 通过自己动手作实验学习如何用插值方法解决实际问题,提高探索和解决问题的能力。通过撰写实验报告,促使自己提炼思想,按逻辑顺序进行整理,并以他人能领会的方式表达自己思想形成的过程和理由。提高写作、文字处理、排版等方面的能力。二、实验内 容 1.编写拉格朗日插值方法的函数M 文件;2.用三种插值方法对已知函数进行插值计算,通过数值和图形输出,比较它们的效果;3.针对实际问题,试建立数学模型,并求解。 三、实验步骤 1.开启软件平台——MATLAB ,开启MATLAB 编辑窗口; 2.根据各种数值解法步骤编写M 文件 3.保存文件并运行; 4.观察运行结果(数值或图形); 5.写出实验报告,并浅谈学习心得体会。 四、实验要求与任务 根据实验内容和步骤,完成以下具体实验,要求写出实验报告(实验目的→问题→数学模型→算法与编程→计算结果→分析、检验和结论→心得体会) 基础实验 1. 一维插值 利用以下一些具体函数,考察分段线性插值、三次样条插值和拉格朗日多项式插值等三种插值方法的差异。 1) 2 11 x +,x ∈[-5,5]; 2)sin x , x ∈[0,2π]; 3)cos 10 x , x ∈[0,2π]. 注意:适当选取节点及插值点的个数;比较时可以采用插值点的函数值与真实函数值的 差异,或采用两个函数之间的某种距离。 2.高维插值 对于二维插值的几种方法:最邻近插值、分片线性插值、双线性插值、三次插值等,利用如下函数进行插值计算,观察其插值效果变化,得出什么结论? 1) ())(sin ),(px t t x f -=ω,参数p =1/2000~1/200;采样步长为:t =4ms~4s ;

计算方法--插值法与拟合实验

实验三 插值法与拟合实验 一、实验目的 1. 通过本实验学会利用程序画出插值函数,并和原图形相比较 2. 通过本实验学会拟合函数图形的画法,并会求平方误差 二、实验题目 1. 插值效果的比较 实验题目:区间[]5,5-10等分,对下列函数分别计算插值节点k x 的值,进行不同类型的插值,作出插值函数的图形并与)(x f y =的图形进行比较: 2 11)(x x f +=; x x f arctan )(=; 4 41)(x x x f += (1) 做拉格朗日插值; (2) 做三次样条插值. 2. 拟合多项式实验 实验题目:给定数据点如下表所示: 分别对上述数据作三次多项式和五次多项式拟合,并求平方误差,作出离散函数),(i i y x 和拟合函数的图形. 三、实验原理 本实验应用了拉格朗日插值程序、三次样条插值程序、多项式拟合程序等实验原理. 四、实验内容 1(1) figure x=-5:0.2:5; y=1./(1+x.^2); plot(x,y,'r'); hold on %拉格朗日插值 x1=-5:1:5; y1=1./(1+x1.^2); xx=-4.5:0.5:4.5; yy=malagr(x1,y1,xx); plot(xx,yy,'+') %三次样条插值 dy0=1./(1+25); dyn=1./(1+25);

m=maspline(x1,y1,dy0,dyn,xx); plot(xx,m,'ok') 1(2) x=-5:0.2:5; y=atan(x); plot(x,y,'r'); hold on %拉格朗日插值 x1=-5:1:5; y1=atan(x1); xx=-4.5:0.5:4.5; yy=malagr(x1,y1,xx); plot(xx,yy,'+') %三次样条插值 dy0=1./(1+25); dyn=1./(1+25); m=maspline(x1,y1,dy0,dyn,xx); plot(xx,m,'ok') 1(3) x=-5:0.2:5; y=x.^2./(1+x.^4); plot(x,y,'r'); hold on %拉格朗日插值 x1=-5:1:5; y1=x1.^2./(1+x1.^4); xx=-4.5:0.5:4.5; yy=malagr(x1,y1,xx); plot(xx,yy,'+') %三次样条插值 dy0=1./(1+25); dyn=1./(1+25); m=maspline(x1,y1,dy0,dyn,xx); plot(xx,m,'ok') 2. x=[-1.5 -1.0 -0.5 0.0 0.5 1.0 1.5]'; y=[-4.45 -0.45 0.55 0.05 -0.44 0.54 4.55]'; plot(x,y,'or'); hold on %三次多项式拟合 p1=mafit(x,y,3);

二元函数插值的一般方法研究

《二元函数多项式插值的一般方法研究》的开题报告 一.课题研究的背景和意义 (一).插值问题的提出和发展过程 许多实际问题都用函数)(x f y =来表示某种内在规律的数量关系,其中相当一部分函数通过实验或观测得到的.虽然)(x f 在某个区间[]b a ,上是存在的,有的还是连续的,但却只能给出[]b a ,上一系列点i x 的函数值),...,1,0)((n i x f y i i ==,这只是一张函数表.有的函数虽有解析表达式,但由于计算复杂,使用不方便,通常也造一个函数表,如大家熟悉的三角函数表、对数表、平方根和立方根表等.为了研究函数的变化规律,往往需要求出不在表上的函数值.因此,我们希望根据给定的函数表做一个既能反应函数)(x f 的特性,又便于计算的简单函数)(x P ,用)(x P 近似)(x f .通常选一类较简单的函数(如代数多项式或分段代数多项式)作为)(x P ,并使)()(i i x f x P =对n i ,...,1,0=成立.这样确定的)(x P 就是我们希望得到的插值函数. 对于上述的)(x f y =的函数插值,前人们已经做过很多的研究,典型的有多项式插值、拉格朗日插值、牛顿插值、埃尔米特插值等.但是对于二元函数),(y x f z =的插值还没有一个较广的研究. (二).二元函数插值研究的意义 1. 理论意义: 一元函数插值主要有基函数法、拉格朗日插值法、牛顿插值法、埃尔米特插值等,但是对于二元函数插值乃至n 元插值是不能直接在一元函数插值的基础上直接推广的。多元插值是一个活跃的研究领域,至今已有非常多的多元插值公式,但是可供利用的公式十分少。 所以我们研究二元函数的插值时,可以为n 元函数插值提供新的研究思路,有助于复杂函数的偏导数的求解,也可以是对插值理论的完善。 2. 实际意义: 一元函数插值问题主要是平面的,而二元函数插值是在三维空间上的,这对我们构造三维空间图像有非常大的作用.例如,在现代机械工业中用计算机控制加工机械零件,根据设

求二元函数极限的几种方法二元函数极限定理

1 / 15 1.二元函数极限概念分析 定义1 设函数f 在2D R ?上有定义,0P 是D 的聚点,A 是一个确定的实数.如果对于任意给定的正数ε,总存在某正数δ,使得00(;)P U P D δ∈时,都有 ()f P A ε-<, 则称f 在D 上当0P P →时,以A 为极限,记0 lim ()P P P D f P A →∈=. 上述极限又称为二重极限. 2.二元函数极限的求法 2.1 利用二元函数的连续性 命题 若函数(,)f x y 在点00(,)x y 处连续,则 0000(,)(,) lim (,)(,)x y x y f x y f x y →=. 例1 求2 (,)2f x y x xy =+ 在点(1,2)的极限. 解: 因为2 (,)2f x y x xy =+在点(1,2)处连续,所以 12 212 2lim (,) lim(2) 12125. x y x y f x y x xy →→→→=+=+??= 例2 求极限()()2 21,1,21 lim y x y x +→. 解: 因函数在()1,1点的邻域内连续,故可直接代入求极限,即 ()()221,1,21lim y x y x +→=31 .

2 / 15 2.2 利用恒等变形法 将二元函数进行恒等变形,例如分母或分子有理化等. 例3 求 00 x y →→ 解 : 00 x y →→ 00 x y →→= 0x y →→= 00 1. 4 x y →→==-例4 ()() 2 2 220,0,321 )31)(21(lim y x y x y x +-++→. 解 : 原式 ()() ( ) )() () ,0,02 211lim 231x y x y →= + ()( 22 ,0,0lim x y →= + 11022 = +=.

相关主题
相关文档 最新文档