%初始化:w10ij表示wij(k-2)；w11ij表示wij(k-1)；w20j表示wjk(k-2)； w21j表示 wjk(k-1)；q和p為閾值
w10ij=[.01 .01 .02; .01 .01 .02; .01 0 .01; .001 .001 .002;.001 0 .002; .0011 .001 .001] ;     
w11ij=[-.1 -.02 .11; -.21 .10 -.19; -.14 .15 -.16; .14 -.13 .17; -.13 .12 .21; -.16 -.23 .13];     
w20j=[.01;.02;.1;.2;.1;.1];  
 w21j=[.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;
k1=1;
p1=.1;
w23j=[0;0;0;0;0;0];w22j=0;
w=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;
n=1; q=0;
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;
k1=1;w0=1;
kn=28;dj2=0.01;
e3=.0055; a1=0.036; a2=0.036; a3=0.08;w0=1;dj2=0.01;
%+++++++++++++++++++++++++++++++++++++初始化結束
for m=1:36
      for n=1:36
 q=pi*0.05*m;
p=pi*0.02*n;
yn=cos(2*q+2*p)*sin(2*p);
 yzs=0.1*yn;
 %xj=[0.4+yzs 0.3 0.2];  %協調器使網絡的輸入隨樣本動態變化
 yp=yn;
    for  k=1:6
        % calculating output of the hidden layer
 for  i=1:6
   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:3
   yb=yb+w21j(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;
      % revising right value 
  for   i=1:6
       d1=[1-o1j(:,i)]*dj2*w23j(i,:);%計算隱層誤差反傳信號
      do=o1j(:,i)*d1;
       q3j(:,i)=q1j(i);
   for j=1:3
      dw=w11ij(i,j)-w10ij(i,j);
    if e2<0.004, a1=0;a2=0;a3=0; else, a1=0.022; a2=0.03;a3=0.003;end;
      w12ij=w11ij(i,j)+a1*do*xj(j)+a2*dw+a3*w;
      w13ij(i,j)=w12ij;
   end
  end      
      w10ij=w11ij;

      w11ij=w13ij;
       q0j=q1j;
      q1j=q3j;
      w=yp-y;  w0=w;
 w=0.36*w0+(yp-y);    %計算積累誤差w  
  if e2<0.004, w=0.78*w; end
  for  i=1:6
      d2j=y*(1-y)*(yp-y);%計算輸出層誤差反傳信號
      dw=w21j(i,:)-w20j(i,:);
      w22j=w21j(i,:)+0.132*d2j*o1j(i)+0.26*dw;+0.026*w;
      w23j(i,:)=w22j;
    end
      w20j=w21j;
      w21j=w23j;
      ph=p1-p0;
      
      p2=p1+.3*(yp-y)+.102*ph;
      p0=p1;
      p1=p2;
     
if  e2<=0.005  break; else  end  % 判斷誤差
end  %對樣本的一次采樣值訓練k次結束
ypp(m,n)=yn;  %存儲二維非線性樣本
yom(m,n)=y;  %存儲神經網絡辨識結果
e3(m,n)=e2;   %存儲訓練誤差
end  %變量n訓練結束
m2=m
end  %變量m訓練結束
w11ij=w13ij  %句末無“；”，可直接在MATLAB主界面下觀察隱層權值
w21j=w23j  %可直接在MATLAB主界面下觀察輸出層權值
%grapher
subplot(2,2,1);
%plot3(m,n,ypp); %也可用此語句繪制樣本的三維圖
mesh(ypp)  %繪制樣本的三維圖
xlabel('m'),ylabel('n'),zlabel('ypp')， %三維坐標
title('Identified model by inp. algorithm')， %圖題
subplot(2,2,2); mesh(yom)
xlabel('m'),ylabel('n'),zlabel('yom');
subplot(2,1,2); mesh(e3)
xlabel('m'),ylabel('n'),zlabel('e3');