算法

1. 雷电

tic;
1x=100;
1y=100;
 a=(cos(pi/lx)+cos(pi/ly))/2; 
w=2/(1+sqrt(1-a*a)); 
v1=zeros(ly,1x); 
v=V1; 
r=3; 
trace=V; 
fori=1:ly
    forj=1:1x
        v1(ij)=(i-1)*100/(1y-1); 
        end
end 
o=sqrt(2);
 v(1,:)=1;
v(ly,:)=1;
V(:,1)=1;
v(:,1x)=1; 
k=0; 
i=2;
j=fix(lx/2); 
while(v(ij)~=1)
    v(i,j)=2;
    trace(ij)=1;
    v1(i,j)=0; 
    V2=V1;
    maxt=1; 
    while(maxt>1e-5)
        k=k+1; 
            maxt=0; 
                for m=2:ly-1
                for n=2:1x-1; 
                    if v(m,n)==0
                        v2(m,n)=v2(m,n)+w*(v1(m,n+1)+v1(m+1,n)+v2(m-1,n)+v2(m,n-1)-4*v2(m,n))/4; 
                        t=abs(v2(m,n)-v1(m,n)); 
                        if(t>maxt)maxt=t; 
                        end
                    end
                end 
        V1=v2; 
    end 
p=0;

for i=2:1y-1
    for j=2:1x-1
        if(v(ij)==2)
            if(v(i-1j)==2)p1=0; else p1=v1(i-1j)"r; end; 
            if(v(i+1j)==2)p2=0; else p2=v1(i+1j)~r; end; 
            if(v(ij-1)==2)p3=0; else p3=v1(ij-1)~r; end; 
            if(v(ij+1)==2)p4=0; else p4=v1(ij+1)~r; end; 
            if(v(i-1j-1)==2)p5=0; else p5=(v1(i-1j-1)/o)~r; end; 
            if(v(i-1j+1)==2)p6=0; else p6=(V1(i-1.j+1)/o)^r; end; 
            if(v(i+1j-1)==2)p7=0; else p7=(v1(i+1.j-1)/o)"r; end; 
            if(v(i+1j+1)==2)p8=0; else p8=(v1(i+1j+1)/o)~r; end; 
            p=p1+p2+p3+p4+p5+p6+p7+p8+p; 
    end 
    end; 
end; 
    g=rand(1,1); 
s=0; 
for i=2:1y-1
    for j=2:1x-1
        if(v(ij)==2)
            if(v(i-1j)==2)  p1=0; else p1=(vl(i-1j)"r/p); end; 
            if(v(i+1j)==2)p2=0; else p2=(v1(i+1j)~r/p); end;
            if(v(ij-1)==2)p3=0; else p3=(v1(ij-1)~r/p); end; 
            if(v(ij+1)==2)p4=0; else p4=(v1(ij+1)~r/p); end; 
            if(v(i-1j-1)==2)p5=0; else p5=((V1(i-1.j-1)/o)~r)/p; end; 
            if(v(i-1j+1)==2)p6=0; else p6=((v1(i-1j+1)/o)~r)/p; end; 
            if(v(i+1j-1)==2)p7=0; else p7=((v1(i+1.j-1)/o)~r)/p; end; 
            if(v(i+1j+1)==2)p8=0; else p8=((V1(i+1j+1)/o)~r)/p; end; 
            s=S+p1+p2+p3+p4+p5+p6+p7+p8; 
            if(q>=s-p1-p2-p3-p4-p5-p6-p7-p8)&&(q<S-p2-p3-p4-p5-p6-p7-p8)
                m=i-1;n=j; 
                end; 
            if(q>=S-p2-p3-p4-p5-p6-p7-p8)&&(q<S-p3-p4-p5-p6-p7-p8)
                m=i+1;n=j;
            end; 
                if(q>=-p3-p4-p5-p6-p7-p8)&&(q<s-p4-p5-p6-p7-p8)
                    m=i;n=j-1; end; 
                if(q>=s-p4-p5-p6-p7-p8)&&(q<=S-p5-p6-p7-p8)
                    m=i;n=j+1; end; 
                if(q>=s-p5-p6-p7-p8)&&(q<=s-p6-p7-p8)
                    m=i-1;n=j-1; end; 
                if(q>=s-p6-p7-p8)&&(q<=s-p7-p8)
                    m=i-1;n=j+1; end; 
                if(q>=S-p7-p8)&&(q<=S-p8)
                    m=i+1;n=j-1; end; 
                if(q>=S-p8)&&(q<=s)
                    m=i+1;n=j+1; end; 
                    end
                end
            end i=m j=n; 
end 
k
toc; 
time=toc; 
im=pcolor(trace)
colormap(cool);
set(im,'LineStyle','none');