stimulate1.m

有限差分

% Fn/n(k-1)=0
% Fn/n(k)= -mun(m)*(Fai(k)-Fai(k+1))/h(m)/(exp(Fai(k)-Fai(k+1))-1)+mun(m-1)*(Fai(k)-Fai(k+1))/h(m-1)/(exp(Fai(k)-Fai(k+1))-1)-(1/2*h(m-1)+1/2*h(m))*p(k)/(tn*(p(k)+1)+tp*(n(k)+1))+(1/2*h(m-1)+1/2*h(m))*(n(k)*p(k)-1)/(tn*(p(k)+1)+tp*(n(k)+1))^2*tp
% Fn/n(k+1)= -(Fai(k)-Fai(k+1))*exp(Fai(k)-Fai(k+1))*(-mun(m)*h(m-1)+mun(m-1)*h(m))/h(m)/(exp(Fai(k)-Fai(k+1))-1)/h(m-1)
% 
% Fn/f(k-1)=0
% Fn/f(k)= -(-mun(m)*h(m-1)+mun(m-1)*h(m))*(Fai(k+1)*n(k+1)*exp(Fai(k)-Fai(k+1))+n(k+1)*exp(Fai(k)-Fai(k+1))^2-n(k)*exp(Fai(k)-Fai(k+1))-n(k+1)*exp(Fai(k)-Fai(k+1))-Fai(k)*n(k+1)*exp(Fai(k)-Fai(k+1))-Fai(k+1)*n(k)*exp(Fai(k)-Fai(k+1))+Fai(k)*n(k)*exp(Fai(k)-Fai(k+1))+n(k))/h(m)/(exp(Fai(k)-Fai(k+1))-1)^2/h(m-1)
%  Fn/f(k+1)= (-mun(m)*h(m-1)+mun(m-1)*h(m))*(n(k+1)*exp(Fai(k)-Fai(k+1))^2+Fai(k+1)*n(k+1)*exp(Fai(k)-Fai(k+1))+Fai(k)*n(k)*exp(Fai(k)-Fai(k+1))-Fai(k+1)*n(k)*exp(Fai(k)-Fai(k+1))-Fai(k)*n(k+1)*exp(Fai(k)-Fai(k+1))-n(k)*exp(Fai(k)-Fai(k+1))-n(k+1)*exp(Fai(k)-Fai(k+1))+n(k))/h(m)/(exp(Fai(k)-Fai(k+1))-1)^2/h(m-1)
% Fn/p(k)= -1/2*(h(m-1)+h(m))*(n(k)+1)*(tp*n(k)+tn)/(tn*p(k)+tn+tp*n(k)+tp)^2
% Fn/p(k+1)=0
% 
% Fp/n(k-1)=0
% Fp/n(k)= 1/2*(h(m-1)+h(m))*(p(k)+1)*(tn*p(k)+tp)/(tn*p(k)+tn+tp*n(k)+tp)^2
% Fp/n(k+1)=0
% Fp/f(k-1)=0
% Fp/f(k)= (-mun(m)*h(m-1)+mun(m-1)*h(m))*(-p(k)*exp(Fai(k)-Fai(k+1))^2+p(k)*exp(Fai(k)-Fai(k+1))+p(k+1)*exp(Fai(k)-Fai(k+1))-p(k+1)+Fai(k)*p(k)*exp(Fai(k)-Fai(k+1))-Fai(k+1)*p(k)*exp(Fai(k)-Fai(k+1))-Fai(k)*p(k+1)*exp(Fai(k)-Fai(k+1))+Fai(k+1)*p(k+1)*exp(Fai(k)-Fai(k+1)))/h(m)/(exp(Fai(k)-Fai(k+1))-1)^2/h(m-1)
% Fp/f(k+1)= -(-mun(m)*h(m-1)+mun(m-1)*h(m))*(-p(k)*exp(Fai(k)-Fai(k+1))^2+p(k)*exp(Fai(k)-Fai(k+1))+p(k+1)*exp(Fai(k)-Fai(k+1))-p(k+1)+Fai(k)*p(k)*exp(Fai(k)-Fai(k+1))-Fai(k+1)*p(k)*exp(Fai(k)-Fai(k+1))-Fai(k)*p(k+1)*exp(Fai(k)-Fai(k+1))+Fai(k+1)*p(k+1)*exp(Fai(k)-Fai(k+1)))/h(m)/(exp(Fai(k)-Fai(k+1))-1)^2/h(m-1)
% 
% Fp/p(k-1)=0
% Fp/p(k)= mun(m)*(Fai(k)-Fai(k+1))*exp(Fai(k)-Fai(k+1))/h(m)/(exp(Fai(k)-Fai(k+1))-1)-mun(m-1)*(Fai(k)-Fai(k+1))*exp(Fai(k)-Fai(k+1))/h(m-1)/(exp(Fai(k)-Fai(k+1))-1)+(1/2*h(m-1)+1/2*h(m))*n(k)/(tn*(p(k)+1)+tp*(n(k)+1))-(1/2*h(m-1)+1/2*h(m))*(n(k)*p(k)-1)/(tn*(p(k)+1)+tp*(n(k)+1))^2*tn
% Fp/p(k+1)= (Fai(k)-Fai(k+1))*(-mun(m)*h(m-1)+mun(m-1)*h(m))/h(m)/(exp(Fai(k)-Fai(k+1))-1)/h(m-1)
%   
% Ff/f(k-1)=1/((2*(h(m-1)+h(m)))*h(m-1))
% Ff/f(k)= -(1/(2*(h(m-1)+h(m)))*h(m-1))*(1/h(m-1)+1/h(m))
% Ff/f(k+1)= 1/((2*(h(m-1)+h(m)))*h(m))
% Ff/n(k)=-1
% Ff/p(k)=1

% Fn =
%  
% mun(m)*(Fai(k)-Fai(k+1))*(n(k+1)*exp(Fai(k)-Fai(k+1))-n(k))/h(m)/(exp(Fai(k)-Fai(k+1))-1)-g*(Fai(k)-Fai(k+1))*(n(k+1)*exp(Fai(k)-Fai(k+1))-n(k))/h(m-1)/(exp(Fai(k)-Fai(k+1))-1)-(1/2*h(m-1)+1/2*h(m))*(n(k)*p(k)-1)/(tn*(p(k)+1)+tp*(n(k)+1))
%  
%  
% Fp =
%  
% mun(m)*(Fai(k)-Fai(k+1))*(p(k)*exp(Fai(k)-Fai(k+1))-p(k+1))/h(m)/(exp(Fai(k)-Fai(k+1))-1)-g*(Fai(k)-Fai(k+1))*(p(k)*exp(Fai(k)-Fai(k+1))-p(k+1))/h(m-1)/(exp(Fai(k)-Fai(k+1))-1)+(1/2*h(m-1)+1/2*h(m))*(n(k)*p(k)-1)/(tn*(p(k)+1)+tp*(n(k)+1))
% R(k) =
%  
% (n(k)*p(k)-1)/(tn*(p(k)+1)+tp*(n(k)+1))
%Ff=
% 1/((2*(h(m-1)+h(m)))*h(m-1))*Fai(k-1)-(1/(2*(h(m-1)+h(m)))*h(m-1))*(1/h(m-1)+1/h(m))*Fai(k)+ 1/((2*(h(m-1)+h(m)))*h(m))*Fai(k+1)+p(k)-n(k)+N(k)
P(k)=[0 0 0;0 0 0;0 0 1/((2*(h(m-1)+h(m)))*h(m-1));]
Q(k)=[-mun(m)*(Fai(k)-Fai(k+1))/h(m)/(exp(Fai(k)-Fai(k+1))-1)+mun(m-1)*(Fai(k)-Fai(k+1))/h(m-1)/(exp(Fai(k)-Fai(k+1))-1)-(1/2*h(m-1)+1/2*h(m))*p(k)/(tn*(p(k)+1)+tp*(n(k)+1))+(1/2*h(m-1)+1/2*h(m))*(n(k)*p(k)-1)/(tn*(p(k)+1)+tp*(n(k)+1))^2*tp, -1/2*(h(m-1)+h(m))*(n(k)+1)*(tp*n(k)+tn)/(tn*p(k)+tn+tp*n(k)+tp)^2, -(-mun(m)*h(m-1)+mun(m-1)*h(m))*(Fai(k+1)*n(k+1)*exp(Fai(k)-Fai(k+1))+n(k+1)*exp(Fai(k)-Fai(k+1))^2-n(k)*exp(Fai(k)-Fai(k+1))-n(k+1)*exp(Fai(k)-Fai(k+1))-Fai(k)*n(k+1)*exp(Fai(k)-Fai(k+1))-Fai(k+1)*n(k)*exp(Fai(k)-Fai(k+1))+Fai(k)*n(k)*exp(Fai(k)-Fai(k+1))+n(k))/h(m)/(exp(Fai(k)-Fai(k+1))-1)^2/h(m-1);
      1/2*(h(m-1)+h(m))*(p(k)+1)*(tn*p(k)+tp)/(tn*p(k)+tn+tp*n(k)+tp)^2,mun(m)*(Fai(k)-Fai(k+1))*exp(Fai(k)-Fai(k+1))/h(m)/(exp(Fai(k)-Fai(k+1))-1)-mun(m-1)*(Fai(k)-Fai(k+1))*exp(Fai(k)-Fai(k+1))/h(m-1)/(exp(Fai(k)-Fai(k+1))-1)+(1/2*h(m-1)+1/2*h(m))*n(k)/(tn*(p(k)+1)+tp*(n(k)+1))-(1/2*h(m-1)+1/2*h(m))*(n(k)*p(k)-1)/(tn*(p(k)+1)+tp*(n(k)+1))^2*tn,(-mun(m)*h(m-1)+mun(m-1)*h(m))*(-p(k)*exp(Fai(k)-Fai(k+1))^2+p(k)*exp(Fai(k)-Fai(k+1))+p(k+1)*exp(Fai(k)-Fai(k+1))-p(k+1)+Fai(k)*p(k)*exp(Fai(k)-Fai(k+1))-Fai(k+1)*p(k)*exp(Fai(k)-Fai(k+1))-Fai(k)*p(k+1)*exp(Fai(k)-Fai(k+1))+Fai(k+1)*p(k+1)*exp(Fai(k)-Fai(k+1)))/h(m)/(exp(Fai(k)-Fai(k+1))-1)^2/h(m-1);
      -1,1,-(1/(2*(h(m-1)+h(m)))*h(m-1))*(1/h(m-1)+1/h(m));]
  W(k)=[-(Fai(k)-Fai(k+1))*exp(Fai(k)-Fai(k+1))*(-mun(m)*h(m-1)+mun(m-1)*h(m))/h(m)/(exp(Fai(k)-Fai(k+1))-1)/h(m-1),0,(-mun(m)*h(m-1)+mun(m-1)*h(m))*(n(k+1)*exp(Fai(k)-Fai(k+1))^2+Fai(k+1)*n(k+1)*exp(Fai(k)-Fai(k+1))+Fai(k)*n(k)*exp(Fai(k)-Fai(k+1))-Fai(k+1)*n(k)*exp(Fai(k)-Fai(k+1))-Fai(k)*n(k+1)*exp(Fai(k)-Fai(k+1))-n(k)*exp(Fai(k)-Fai(k+1))-n(k+1)*exp(Fai(k)-Fai(k+1))+n(k))/h(m)/(exp(Fai(k)-Fai(k+1))-1)^2/h(m-1);
       0,(Fai(k)-Fai(k+1))*(-mun(m)*h(m-1)+mun(m-1)*h(m))/h(m)/(exp(Fai(k)-Fai(k+1))-1)/h(m-1),-(-mun(m)*h(m-1)+mun(m-1)*h(m))*(-p(k)*exp(Fai(k)-Fai(k+1))^2+p(k)*exp(Fai(k)-Fai(k+1))+p(k+1)*exp(Fai(k)-Fai(k+1))-p(k+1)+Fai(k)*p(k)*exp(Fai(k)-Fai(k+1))-Fai(k+1)*p(k)*exp(Fai(k)-Fai(k+1))-Fai(k)*p(k+1)*exp(Fai(k)-Fai(k+1))+Fai(k+1)*p(k+1)*exp(Fai(k)-Fai(k+1)))/h(m)/(exp(Fai(k)-Fai(k+1))-1)^2/h(m-1);
       0,0,1/((2*(h(m-1)+h(m)))*h(m));]
   H(k)=[-( mun(m)*(Fai(k)-Fai(k+1))*(n(k+1)*exp(Fai(k)-Fai(k+1))-n(k))/h(m)/(exp(Fai(k)-Fai(k+1))-1)-g*(Fai(k)-Fai(k+1))*(n(k+1)*exp(Fai(k)-Fai(k+1))-n(k))/h(m-1)/(exp(Fai(k)-Fai(k+1))-1)-(1/2*h(m-1)+1/2*h(m))*(n(k)*p(k)-1)/(tn*(p(k)+1)+tp*(n(k)+1)));
         -(mun(m)*(Fai(k)-Fai(k+1))*(p(k)*exp(Fai(k)-Fai(k+1))-p(k+1))/h(m)/(exp(Fai(k)-Fai(k+1))-1)-g*(Fai(k)-Fai(k+1))*(p(k)*exp(Fai(k)-Fai(k+1))-p(k+1))/h(m-1)/(exp(Fai(k)-Fai(k+1))-1)+(1/2*h(m-1)+1/2*h(m))*(n(k)*p(k)-1)/(tn*(p(k)+1)+tp*(n(k)+1)));
         -( 1/((2*(h(m-1)+h(m)))*h(m-1))*Fai(k-1)-(1/(2*(h(m-1)+h(m)))*h(m-1))*(1/h(m-1)+1/h(m))*Fai(k)+ 1/((2*(h(m-1)+h(m)))*h(m))*Fai(k+1)+p(k)-n(k)+N(k));]