双曲偏微分方程算法.txt

求解偏微分方程的matlab程序

%%双曲型偏微分方程组数值解法  例：Ut+AUx=0
A1=-[0.4999840,0.1718901,0.1055926,0.0363018
    0.1718901,0.4166958,0.0363018,0.0880028
    0.1055926,0.0363018,0.4425978,0.1521612
    0.0363018,0.0880028,0.1521612,0.3688691];
B=-[0.4999836,-0.1038723,0.1055925,-0.0219370
    -0.1038723,0.5736108,-0.0219370,0.1211420
    0.1055925,-0.0219370,0.4425974,-0.0919502
    -0.0219370,0.1211420,-0.0919502,0.5077740];
A=inv(A1)*B;
[V,D]=eig(A');%%求A的转置的特征值以及A的转置的右特征向量V1
S=V';%A的转置的右特征向量就是A的左特征向量

%以下构造特征值r(i)组成的对角阵
for i=1:4
    r1(i)=max(D(i,i),0);
    r2(i)=-max(D(i,i),0);
end
Azheng=diag(r1);
Afu=diag(r2);
A1=inv(S)*Azheng*S;
A2=inv(S)*Afu*S;
E=eye(4);

%赋初值
t=0.01;h=0.2;
nt=1/t;nh=1/h;r=t/h;
for j=1:(nh+1)
    u(1,j,:)=[0.5;0.5;0.5;0.5];
end
for i=2:(nt+1)
    u(i,1,:)=[0;0;0;0];
    u(i,nh+1,:)=[0;0;0;0];
end

%%%迭代
%w=input('Please input the N: ');
for i=1:nt
    for j=2:nh
        u(i+1,j,:)=reshape(r*A2*reshape(u(i,j+1,:),4,1)+(E-r*(A1+A2))*reshape(u(i,j,:),4,1)+r*A1*reshape(u(i,j-1,:),4,1),4,1);
    end
end
%%检验结果的收敛性
for i=3:10
jieguo=reshape(u(i,3,:),1,4);
disp(jieguo)
end