matlab数值分析程序--高等数学,数值代数的matlab实现-文字版, matlab电子书, 和matlab 有关的电子书:

6.7.4 用一阶导数计算的几种样条函数和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,由 的坐标组成的向量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

运行后输出的值如下

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

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)= 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,由 的坐标组成的向量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   

   285/43   

输入程序

>> x2=3,

s2=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

x4=8,

s4=s4=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

运行后输出的值如下

x2 =    s2 =                     x 4=    s4 =

  3      24.6221         8       68.5581

例6.7.3和例6.7.2用的方法不同,所以计算的结果不同,但是两种方法计算的结果相近.

 

欢迎转载,转载请注明来自一手册:http://yishouce.com/book/3/3060704.html
友情链接It题库(ittiku.com)| 版权归yishouce.com所有| 友链等可联系 admin#yishouce.com|粤ICP备16001685号-1