dya.m

介绍了DYa算法的具体应用,是解决偏微分方程数值解的交换隐格式的方法

function DYa(f,g1,g2,g3,g4,g5,Za,Zb,h1,h2)
tic;
r=h2/h1^2;
m=(Za(2)-Za(1))/h1;
n=(Zb(2)-Zb(1))/h2;
x=Za(1):h1:Za(2);
y=Za(1):h1:Za(2);
z=Zb(1):h2:Zb(2);
U=zeros(m+1,m+1,n+1);
V=zeros(m+1,m+1);
for t=2:m
    U(t,2:m,1)=feval(g1,x(2:m),y(t));
end
for t=1:n+1
    U(1:m+1,1,t)=feval(g2,y,z(t));
    U(1:m+1,m+1,t)=feval(g3,y,z(t));
end
for t=1:n+1
    U(1,2:m,t)=feval(g4,x(2:m),z(t));
    U(m+1,2:m,t)=feval(g5,x(2:m),z(t));
end
A=ones(1,m-1)*(1+r);
C=ones(1,m-2)*(-r/2);
D=ones(1,m-2)*(-r/2);
A1=ones(1,m-1)*(r/2)*(1-r);
C1=ones(1,m-2)*(r^2/4);
D1=ones(1,m-2)*(r^2/4);
A2=ones(1,m-1)*(1-r)^2;
C2=ones(1,m-2)*(r/2)*(1-r);
D2=ones(1,m-2)*(r/2)*(1-r);
for k=1:n
    for t=2:m
        for s=2:m
            if s==2
                B(1)=A1(1)*(U(t-1,2,k)+U(t+1,2,k))+D1(1)*(U(t-1,3,k)+U(t+1,3,k))+A2(1)*U(t,2,k)+D2(1)*U(t,3,k)+h2/2*(feval(f,x(2),y(t),z(k))+feval(f,x(2),y(t),z(k+1)));
            else if s==m
                    B(m-1)=C1(m-2)*(U(t-1,m-1,k)+U(t+1,m-1,k))+A1(m-1)*(U(t-1,m,k)+U(t+1,m,k))+C2(m-2)*U(t,m-1,k)+A2(m-1)*U(t,m,k)+h2/2*(feval(f,x(m),y(t),z(k))+feval(f,x(m),y(t),z(k+1)));
                else
                    B(s-1)=C1(s-2)*(U(t-1,s-1,k)+U(t+1,s-1,k))+A1(s-1)*(U(t-1,s,k)+U(t+1,s,k))+D1(s-1)*(U(t-1,s+1,k)+U(t+1,s+1,k))+C2(s-2)*U(t,s-1,k)+A2(s-1)*U(t,s,k)+D2(s-1)*U(t,s+1,k)+h2/2*(feval(f,x(s),y(t),z(k+1))+feval(f,x(s),y(t),z(k)));
                end
            end
        end
        B(1)=B(1)+r^2/4*(U(t-1,1,k)+U(t+1,1,k))+r/2*(1-r)*U(t,1,k)+r/2*(U(t,1,k+1)-r/2*(U(t-1,1,k+1)-2*U(t,1,k+1)+U(t+1,1,k+1)));
        B(m-1)=B(m-1)+r^2/4*(U(t-1,m+1,k)+U(t+1,m+1,k))+r/2*(1-r)*U(t,m+1,k)+r/2*(U(t,m+1,k+1)-r/2*(U(t-1,m+1,k+1)-2*U(t,m+1,k+1)+U(t+1,m+1,k+1)));
        V(t,2:m)=zhuiganfa(A,C,D,B);
    end
    A3=ones(1,m-1)*(1+r);
    C3=ones(1,m-2)*(-r/2);
    D3=ones(1,m-2)*(-r/2);
    for t=2:m
        B3(1:m-1)=V(2:m,t);
        B3(1)=B3(1)+r/2*U(1,t,k+1);
        B3(m-1)=B3(m-1)+r/2*U(m+1,t,k+1);
        U(2:m,t,k+1)=zhuiganfa(A3,C3,D3,B3);
    end
end
[U(m/4+1,m/4+1,n/4+1),U(3*m/4+1,m/4+1,n/4+1),U(m/4+1,3*m/4+1,n/4+1),U(3*m/4+1,3*m/4+1,n+1)]
toc