www.pudn.com > MATLABpro.rar > xycf.m

%y''=p(x)y'+q(x)y+r(x),a<=x<=b 
%y(a)=af,y(b)=bt; 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%实例%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
  p=inline('-2/x'); 
  q=inline('2/x^2'); 
  r=inline('sin(log10(x)/log10(exp(1)))/x^2'); 
  N=9; 
  a0=1;b0=2; 
  af=1;bt=2; 
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
  h=(b0-a0)/(N+1); 
  x=a0+h; 
  a(1)=2+h*h*q(x); 
  b(1)=-1+(h/2)*p(x); 
  d(1)=-h*h*r(x)+(1+(h/2)*p(x))*af; 
    for i=2:N-1 
        x=a0+i*h; 
        a(i)=2+h*h*q(x); 
        b(i)=-1+(h/2)*p(x); 
        c(i)=-1-(h/2)*p(x); 
        d(i)=-h*h*r(x); 
     end 
     x=b0-h; 
     a(N)=2+h*h*q(x); 
     c(N)=-1-(h/2)*p(x); 
     d(N)=-h*h*r(x)+(1-(h/2)*p(x))*bt; 
     %%%%%%%%%追赶法%%%%%%%%%%%%%%%%%% 
     %y=trisys(c,a,b,d) 
     L(1)=a(1); 
     u(1)=b(1)/a(1); 
     for i=2:N-1 
         L(i)=a(i)-c(i)*u(i-1); 
         u(i)=b(i)/L(i); 
     end 
     L(N)=a(N)-c(N)*u(N-1); 
     z(1)=d(1)/L(1); 
     for i=2:N 
         z(i)=(d(i)-c(i)*z(i-1))/L(i); 
     end 
     y(N)=z(N); 
     for i=N-1:-1:1 
         y(i)=z(i)-u(i)*y(i+1); 
     end 
     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
     Y=[af,y,bt]; 
 for i=1:N+2 
      x=a0+(i-1)*h; 
      zj(i)=1.1392070132*x-0.03920701320/x^2-3*sin(log10(x)/log10(exp(1)))/10-cos(log10(x)/log10(exp(1)))/10; 
 end 
disp('下面两列分别是数值解和近似解'); 
 re=[Y'      zj']