[½»Á÷] 䶨Òå±äÁ¿W£¬ÇóÖú£¡£¡£¡Ð»Ð» ÒÑÓÐ2È˲ÎÓë

function DP=f2(t,P)
global A N
format long g;
N=15;
S=zeros(1,N);
[m,l]=gausspoint(N-2);
S(1)=-1;
S(N)=1;
for i=2:N-1
    S(i)=m(i-1);
end
S;
X=zeros(1,N);
X(1)=0;
X(2)=3*10^(-4);
X(N-1)=1-5*10^(-5);
X(N)=1;
for i=3:N-2
   X(i)=0.5*S(i)+0.5;
end
%ÇóȨϵÊý
A=zeros(N,N,4);
XX=ones(1,N);
for i=1:N
    for j=1:i-1
        XX(i)=(X(i)-X(j))*XX(i);
    end
    for j=i+1:N
        XX(i)=(X(i)-X(j))*XX(i);
    end
end
for i=1:N
    for j=1:i-1
        A(i,j,1)=XX(i)/XX(j)/(X(i)-X(j));
    end
    for j=i+1:N
        A(i,j,1)=XX(i)/XX(j)/(X(i)-X(j));
    end
end
QM=-sum(A(:,:,1)');
for i=1:N
    A(i,i,1)=QM(i);
end
for r=2:4
    for i=1:N
        for j=1:i-1
            A(i,j,r)=r*(A(i,i,r-1)*A(i,j,1)-A(i,j,r-1)/(X(i)-X(j)));
        end
        for j=i+1:N            
            A(i,j,r)=r*(A(i,i,r-1)*A(i,j,1)-A(i,j,r-1)/(X(i)-X(j)));
        end
    end
   HM= -sum(A(:,:,r)');
   for i=1:N
       A(i,i,r)=HM(i);
   end
end
r1=0;r2=40;v=0.25
for i=1:N
    for k=1:N
        if i==2
        K(i,k)=A(i,k,2)+v/(X(i))*A(i,k,1);
        elseif i==N-1;
        K(i,k)=A(i,k,2)+v/(X(i))*A(i,k,1);
        else
          K(i,k)=A(i,k,4)+2/(X(i))*A(i,k,2)-1/(X(i)^2)*A(i,k,2)+1/(X(i)^3)*A(i,k,1);  
        end
    end
end
K1=-K(2:N-1,2:N-1);%B¾ØÕó
K2=eye(2*N-2);
K=zeros(26,26);
for i=1:N-2
    for j=1:N-2
  K=K1;
    end;
end;
for i=N-1:(2*N-4)
    for j=N-1:(2*N-4)
        if i==j
  K(i,j)=1;
        end;
    end;
end;
DL=K1*W;
DW=K2*L;
Q=[W;L];
P=[L;W];
DP=K*Q;
end

2Â¥: Originally posted by Ô¼ÊøÌÕ´É×°¼× at 2017-04-17 11:23:08
WÊÇɶ£¿LÊÇɶ£¿

3Â¥: Originally posted by quarkhe at 2017-04-17 18:02:33
ÊÇÁ½¸öδ֪´ýÇó½âµÄÁоØÕó...

3Â¥: Originally posted by quarkhe at 2017-04-17 18:02:33
ÊÇÁ½¸öδ֪´ýÇó½âµÄÁоØÕó...

4Â¥: Originally posted by Ô¼ÊøÌÕ´É×°¼× at 2017-04-17 19:35:39
¼ÈÈ»ÊÇδ֪Á¿£¬ÔõôÄܷŵ½¸³ÖµºÅÓұߣ¿...