clear
w10ij=[.001 .001 .002; .001 .001 .02; .01 0.01 .001; .001 .001 .002;.001 0.01 .002; .0011 .001 .001]     
w11ij=[-.1 -.02 .11; -.21 .10 -.19; -.14 .15 -.16; .14 -.13 .17; -.13 .12 .21; -.16 -.23 .13]    
w20jk=[.01;.02;.1;.2;.1;.1] 
 w21jk=[.18;.9;.9;.7;.8;.9]   
q0j=[.5 .8 .4 .6 .1 .2 ]    
q1j=[-.1 .02 .12 .14 -.02 .02];   
q2j=0;                      % threshold value
p0=.2;
p1=.2;
w23jk=[0;0;0;0;0;0];w22jk=0;
 xj=[0.5 0.3 0.2];           % inputs 
ya=[0 0 0]; yp=0; yy=0; m1=0;yam=0;yp1=0; qw=0; 
yo=[0 0 0]; ya1=0;
error=0.0001;
e1=0; e0=0;  e2=0;  e3=0; e4=0; 
yo=0; ya=0; yb=0; y0=0; y1=0; y2=0; y3=0; 
u=0; u1=0; u2=0;
dj2=0
e3=0; a1=0.02; a2=0.05; a3=0.005; M=120;  dj2=0;
%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
for m=1:M    
p=pi*0.02*m;
yn=sin(2*p);
 yp=yn;
    for  k=1:4
 for  i=1:6    % calculating output of the hidden layer
   x1=[w11ij(i,1)*xj(:,1)]+[w11ij(i,2)*xj(:,2)]+[w11ij(i,3)*xj(:,3)];
   x=x1+q1j(:,i);
   o=tanh(x);
   o1j(i)=o;
  end  
 % calculating output of the output layer
   for  i=1:6
   yb=yb+w21jk(i,:)*o1j(:,i);
   end
   yi=yb+p1;
   y=tanh(yi);
 % calculating error value between aim and practice value  
   e0=e1;
   e1=e2;
   e2=[(yp-y).^2]/2;
 
  for   i=1:6     
        d1=[1-o1j(:,i)]*dj2*w23jk(i,:); %計算隱層誤差反傳信號
      do=o1j(:,i)*d1;
       qw=q1j(:,i)-q0j(:,i);
      q2j=q1j(:,i)+.12*do+.01*qw;
       q3j(:,i)=q2j;
   for j=1:3
      dw=w11ij(i,j)-w10ij(i,j);
      if e2<0.05, a1=0;a2=0;a3=0; else, a1=0.032; a2=0.05;a3=0.005;end;
       w12ij=w11ij(i,j)+a1*do*xj(:,j)+a2*dw;
       
      w13ij(i,j)=w12ij;
   end
  end      
      w10ij=w11ij;
      w11ij=w13ij;
       q0j=q1j;
      q1j=q3j;
  for  i=1:6
      d2j=y*(1-y)*(yp-y);%計算輸出層誤差反傳信號
      dw=w21jk(i,:)-w20jk(i,:);
      w22jk=w21jk(i,:)+0.1*d2j*o1j(i)+0.11*dw;
      w23jk(i,:)=w22jk;
    end
      w20jk=w21jk;
      w21jk=w23jk;
      ph=p1-p0;
      p2=p1+.6*(yp-y)+.08*ph;
      p0=p1;
      p1=p2;
if  e2<=.005  break;
else
end
end
ypp(m)=yn;
yom(m)=y;
e3(m)=e2;
end, 
w13ij   %打印隱層權值；
w23jk   %打印輸出層權值；
%grapher
m=1:M;
subplot(4,1,1);
plot(ypp)
%xlabel('m'), %若去掉該行最前的“%”，第一幅圖的橫坐標軸會標出“m”；
ylabel('ypp'),
%title('Identified model by BP. algorithm')
subplot(4,1,2);
plot(m,yom,m,yom,'rx')
%xlabel('m'),
ylabel('yom'),
subplot(4,1,3);
plot(e3)
%xlabel('m'),
ylabel('e3'),
subplot(4,1,4);
plot(m,ypp,m,yom,'gx')
xlabel('m'),ylabel('ypp&yom'),
legend('ypp=sin(2*pi*m/50)','yom=output of NN'); %圖標注