主题:[转帖]MATLAB解偏微分方程
如有问题联系:1759553924
function [err,zuidaerr,U,Ujing]=houeulererwei(m,n,jingqueu,phi,f)
%此程序用于求解差分方程du/dt-d2u(x)-d2u(y)=f(x)的解
%u(x,0)=phi,u(a,t)=alpha,u(b,t)=beta
%差分空间剖分m等分,时间剖分n等分
h=(1-0)/m;
t=(1-0)/n;
row=t/h;
r=t/h.^2;
x=0:h:1;
alpha=0.001*ones(1,m+1);
y=0:t:1;
for i=1:m+1
U(i,1:m+1)=feval(phi,(i-1)*h,x);
end
U(1:m+1,1)=zeros(m+1,1);
U(1:m+1,1+m)=zeros(m+1,1);
U(1,1:m+1)=zeros(1,1+m);
U(m+1,1:1+m)=zeros(1,1+m);
A=(1/h/2)*ones(1,m-2)-[alpha(2:m-1)+alpha(3:m)]/2/h^2;
C=-1/h/2*ones(1,m-2)-[alpha(2:m-1)+alpha(3:m)]/2/h^2;
B=(1/t)*ones(1,m-1)+[2*alpha(2:m)+alpha(1:m-1)+alpha(3:m+1)]/h^2;
AA=diag(B)+diag(C,-1)+diag(A,1);
AAAA=[];
for k=1:m-1
AAAA=blkdiag(AAAA,AA);
end
A1=(1/h/2)*ones(1,m-1)-[alpha(2:m)+alpha(3:m+1)]/2/h^2;
BB=diag(A1);
BBBB=[];
for k=1:m-2
BBBB=blkdiag(BBBB,BB);
end
BBBB=[BBBB;zeros(m-1,(m-1)*(m-2))];
BBBB=[zeros((m-1)^2,m-1) BBBB];
C2=-1/h/2*ones(1,m-1)-[alpha(1:m-1)+alpha(2:m)]/2/h^2;
CC=diag(C2);
CCCC=[];
for k=1:m-2
CCCC=blkdiag(CCCC,CC);
end
CCCC=[CCCC zeros((m-1)*(m-2),m-1)];
CCCC=[zeros(m-1,(m-1)^2) ; CCCC];
DDDD=AAAA+BBBB+CCCC;
for k=1:n
ff=[];
for j=1:m-1
ff=[ff 1/t*U(j+1,2:m)+[feval(f,h:h:1-h,h*j,k*t)]];
end
xxx=DDDD\ff';
xxx=reshape(xxx,m-1,m-1);
U(2:m,2:m)=xxx';
for j=1:m+1
Ujing(j,1:m+1)=feval(jingqueu,(j-1)*h,x,t*k);
end
zuidaerr(k)=max(max(abs(U-Ujing)))
end
err=abs(U-Ujing);
%显示四位有效数
%zuidaERR=vpa(zuidaERR,4);
%k=n/10+1:n/10:n+1;
%U=[U(k,6)];
%U=vpa(U,7);
%ERR=[ERR(k,6)];
%ERR=vpa(ERR,4);
function [err,zuidaerr,U,Ujing]=houeulererwei(m,n,jingqueu,phi,f)
%此程序用于求解差分方程du/dt-d2u(x)-d2u(y)=f(x)的解
%u(x,0)=phi,u(a,t)=alpha,u(b,t)=beta
%差分空间剖分m等分,时间剖分n等分
h=(1-0)/m;
t=(1-0)/n;
row=t/h;
r=t/h.^2;
x=0:h:1;
alpha=0.001*ones(1,m+1);
y=0:t:1;
for i=1:m+1
U(i,1:m+1)=feval(phi,(i-1)*h,x);
end
U(1:m+1,1)=zeros(m+1,1);
U(1:m+1,1+m)=zeros(m+1,1);
U(1,1:m+1)=zeros(1,1+m);
U(m+1,1:1+m)=zeros(1,1+m);
A=(1/h/2)*ones(1,m-2)-[alpha(2:m-1)+alpha(3:m)]/2/h^2;
C=-1/h/2*ones(1,m-2)-[alpha(2:m-1)+alpha(3:m)]/2/h^2;
B=(1/t)*ones(1,m-1)+[2*alpha(2:m)+alpha(1:m-1)+alpha(3:m+1)]/h^2;
AA=diag(B)+diag(C,-1)+diag(A,1);
AAAA=[];
for k=1:m-1
AAAA=blkdiag(AAAA,AA);
end
A1=(1/h/2)*ones(1,m-1)-[alpha(2:m)+alpha(3:m+1)]/2/h^2;
BB=diag(A1);
BBBB=[];
for k=1:m-2
BBBB=blkdiag(BBBB,BB);
end
BBBB=[BBBB;zeros(m-1,(m-1)*(m-2))];
BBBB=[zeros((m-1)^2,m-1) BBBB];
C2=-1/h/2*ones(1,m-1)-[alpha(1:m-1)+alpha(2:m)]/2/h^2;
CC=diag(C2);
CCCC=[];
for k=1:m-2
CCCC=blkdiag(CCCC,CC);
end
CCCC=[CCCC zeros((m-1)*(m-2),m-1)];
CCCC=[zeros(m-1,(m-1)^2) ; CCCC];
DDDD=AAAA+BBBB+CCCC;
for k=1:n
ff=[];
for j=1:m-1
ff=[ff 1/t*U(j+1,2:m)+[feval(f,h:h:1-h,h*j,k*t)]];
end
xxx=DDDD\ff';
xxx=reshape(xxx,m-1,m-1);
U(2:m,2:m)=xxx';
for j=1:m+1
Ujing(j,1:m+1)=feval(jingqueu,(j-1)*h,x,t*k);
end
zuidaerr(k)=max(max(abs(U-Ujing)))
end
err=abs(U-Ujing);
%显示四位有效数
%zuidaERR=vpa(zuidaERR,4);
%k=n/10+1:n/10:n+1;
%U=[U(k,6)];
%U=vpa(U,7);
%ERR=[ERR(k,6)];
%ERR=vpa(ERR,4);