matlab 差分方程画图,Matlab差分方程画图的问题

方程大概是这样的:

x_i (k+1)=f(x_i (k))+g(x_i (k-d(k)))+sum_{m=1}^{infinity}(u_m  *h(x_i (k-m)))+sum_{j=1,j~=i}^{N}w_{ij} x_j (k)

是一个耦合的网络,不妨取N=3,即是3个网络耦合.下面是我编的一个画图程序,大家帮助看一下哪里错了,因为总是和手算的值不一样,这里的m的上界本应该是无穷,但是那样比较麻烦,我只取到1000,应该和实际比较接近才对.

function disc_nn()

tic; % 开始计时

clear;

k0=0; kfinal=50; h=1;

i=1;

for k=k0:h:kfinal % time range

d(i)=3+(1+(-1)^k)/2;

i=i+1;

end

u1(1,:)=fun1_ini(0);

u2(1,:)=fun2_ini(0);

u3(1,:)=fun3_ini(0);

i=2;

sumh11=0;

sumh12=0;

sumh21=0;

sumh22=0;

sumh31=0;

sumh32=0;

for k=k0+1:h:kfinal % time range

for m=1:h:1000

if (k-m)<=0 % 判断k-m的符号

u11(k,1)=u1(1,1);

u12(k,2)=u1(1,2);

u21(k,1)=u2(1,1);

u22(k,2)=u2(1,2);

u31(k,1)=u3(1,1);

u32(k,2)=u3(1,2);

else

u11(k,1)=u1((k-m),1);

u12(k,2)=u1((k-m),2);

u21(k,1)=u2((k-m),1);

u22(k,2)=u2((k-m),2);

u31(k,1)=u3((k-m),1);

u32(k,2)=u3((k-m),2);

end

sumh11=sumh11+(2^(-m-3))*h1(u11(k,1));

sumh12=sumh12+(2^(-m-3))*h2(u12(k,2));

sumh21=sumh21+(2^(-m-3))*h1(u21(k,1));

sumh22=sumh22+(2^(-m-3))*h2(u22(k,2));

sumh31=sumh31+(2^(-m-3))*h1(u31(k,1));

sumh32=sumh32+(2^(-m-3))*h2(u32(k,2));

end

if floor(k-d(k))<=0 % 判断k-d(k)近似整数的符号

u1(i,1)=0.05*u2(k,1)+0.05*u3(k,1)...

+f1(u1(k,1))+0.2*u1(k,2)+g1(u1(1,1))+sumh11;

u1(i,2)=0.05*u2(k,2)+0.05*u3(k,2)...

+f2(u1(k,2))+g2(u1(1,2))+sumh12;

u2(i,1)=0.05*u1(k,1)+0.05*u3(k,1)...

+f1(u2(k,1))+0.2*u2(k,2)+g1(u2(1,1))+sumh21;

u2(i,2)=0.05*u1(k,2)+0.05*u3(k,2)...

+f2(u2(k,2))+g2(u2(1,2))+sumh22;

u3(i,1)=0.05*u1(k,1)+0.05*u2(k,1)...

+f1(u3(k,1))+0.2*u3(k,2)+g1(u3(1,1))+sumh31;

u3(i,2)=0.05*u1(k,2)+0.05*u2(k,2)...

+f2(u3(k,2))+g2(u3(1,2))+sumh32;

else

u1(i,1)=0.05*u2(k,1)+0.05*u3(k,1)...

+f1(u1(k,1))+0.2*u1(k,2)+g1(u1(floor(k-d(k)),1))+sumh11;

u1(i,2)=0.05*u2(k,2)+0.05*u3(k,2)...

+f2(u1(k,2))+g2(u1(floor(k-d(k)),2))+sumh12;

u2(i,1)=0.05*u1(k,1)+0.05*u3(k,1)...

+f1(u2(k,1))+0.2*u2(k,2)+g1(u2(floor(k-d(k)),1))+sumh21;

u2(i,2)=0.05*u1(k,2)+0.05*u3(k,2)...

+f2(u2(k,2))+g2(u2(floor(k-d(k)),2))+sumh22;

u3(i,1)=0.05*u1(k,1)+0.05*u2(k,1)...

+f1(u3(k,1))+0.2*u3(k,2)+g1(u3(floor(k-d(k)),1))+sumh31;

u3(i,2)=0.05*u1(k,2)+0.05*u2(k,2)...

+f2(u3(k,2))+g2(u3(floor(k-d(k)),2))+sumh32;

end

i=i+1;

end

%%%%%%%%%%%%%%%%%%%%%%% plot %%%%%%%%%%%%%%%%%%%%%%%%%

n=[-1:h:kfinal];

x11=[55;u1(:,1)];

x12=[-80;u1(:,2)];

x21=[20;u2(:,1)];

x22=[145;u2(:,2)];

x31=[52;u3(:,1)];

x32=[40;u3(:,2)];

figure;

plot(n,x11,'r:.');ylabel('x1');%axis([-1 kfinal -30 60]);

hold on;

plot(n,x21,'b:.');

hold off;

figure;

plot(n,x12,'r:.');ylabel('x2');%axis([-1 kfinal -100 250]);

hold on;

plot(n,x22,'b:.');

hold off;

figure;

plot(n,x11,'r:.');ylabel('x1');%axis([-1 kfinal -30 60]);

hold on;

plot(n,x31,'b:.');

hold off;

figure;

plot(n,x12,'r:.');ylabel('x2');%axis([-1 kfinal -100 250]);

hold on;

plot(n,x32,'b:.');

hold off;

% u1(40:50,:)

% u2(40:50,:)

% u3(40:50,:)

elapse_time=toc % 返回程序的运行时间

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function xini_1=fun1_ini(s)

xini_1(1)=55;

xini_1(2)=-80;

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function xini_2=fun2_ini(s)

xini_2(1)=20;

xini_2(2)=145;

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function xini_3=fun3_ini(s)

xini_3(1)=52;

xini_3(2)=40;

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function f1=f1(x)

f1=tanh(0.2*x)-0.5*x;

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function f2=f2(x)

f2=0.95*x-tanh(0.75*x);

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function g1=g1(x)

g1=0.2*x-tanh(0.1*x);

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function g2=g2(x)

g2=0.1*x;

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function h1=h1(x)

h1=0.2*x-tanh(0.1*x);

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function h2=h2(x)

h2=0.1*x;

return