xsolver.m

别人的边界元程序

function zx=xsolver(equ,k,m,ifp,ep)

%% using "ni bo qian" method to calculate the equations 
global xx
global index
global b
%%  k--source point  m--total number of nodes 
k1=k-1;
mk=m-k;
%% if k~=1, the eqution must be reduing the k-th row 
if k==1
else
    for j=1:k1                 %% exchange the main element
        ll=index(j);        
        if ll==j
        else
            c=equ(j);
            equ(j)=equ(ll);
            equ(ll)=c;
        end
    end
%%  reduing the k-th row   
    for i=1:k1
        kx=(mk+1)*(i-1)+1;                        %% there maybe some errors 
        f=equ(i);
        if f==0
        else                                %%if not zero, reducing the row
            for j=k:m
                ix=kx+j-k;
                equ(j)=equ(j)-xx(ix)*f;
            end
            b(k)=b(k)-b(i)*f;
        end
    end
end
%%   select the main element for k-th row
kp=0;
c=ep;
for km=k:m
    f=abs(equ(km));
    if f<=c
    else
        kp=km;
        c=f;
    end
end
%%  error warnings 
if kp~=0 
else
    disp('***Warning, there is singular in the equation!****');
    return;
end
%%
p=1/equ(kp);
b(k)=b(k)*p;
%%
if k==m
    f=b(m);                  %%beginning solving the equations
    for j=1:k1
        i=m-j;
        ll=index(i);
        c=b(ll);
        b(ll)=b(i)-f*xx(i);
        if i==ll
        else
            b(i)=c;
        end
    end
    return
else
%%    normalized by the k-th element
    equ(kp)=equ(k);
    for j=1:mk
        jk=j+k;
        equ(jk)=equ(jk)*p;
    end
%%    reducing the column
    if k==1
    else
        for i=1:k1
            kx=(mk+1)*(i-1)+1;              %%maybe there are some errors here  !!!!!!!
            kpx=kx+kp-k;                %%exchange the element
            f=xx(kpx);
            xx(kpx)=xx(kx);
            for j=1:mk
                jk=j+k;
                jx=kx+j;
                lx=jx-i;
                xx(lx)=xx(jx)-equ(jk)*f;
            end
            b(i)=b(i)-b(k)*f;
        end
    end
%%    restore the new equation to the matrix
    for j=1:mk
        jk=j+k;
        lx=k1*mk+j;
        xx(lx)=equ(jk);
        index(k)=kp;
    end
    return
%%
end