clear all
close all
clc
%example-3,求解温度场；
a=3;b=2;h=0.05;
x=0:h:a;
y=0:h:b;
[X,Y]=meshgrid(x,y);
Z=0.5*(1/sinh(2*pi*b/a)*sinh(2*pi*Y/a).*sin(2*pi*X/a)+...
    1/sinh(4*pi*b/a)*sinh(4*pi*Y/a).*sin(4*pi*X/a))+...
    1/sinh(3*pi*a/b)*sinh(3*pi*X/b).*sin(3*pi*Y/b);
figure;surf(X,Y,Z); 
title('解析法得到的温度场的三维图')
figure;contour(X,Y,Z,20); %控制等高线的条数
title('解析法得到的温度场的等高图')

figure; meshc(X,Y,Z);
colormap hsv
title('解析法得到的温度场的三维图与等高图')
%============数值解-有限差分法===========================
XL=0;XR=a;
YL=0;YR=b;
N=(XR-XL)/h;
M=(YR-YL)/h;
u=zeros(N+1,M+1);
u(N+1,1:M+1)=sin(3*pi/b*(0:M)*h); %边界条件
u(1:N+1,M+1)=sin(3*pi/a*(0:N)*h).*cos(pi/a*(0:N)*h);
uold=u;
u(2:N,2:M)=u(2:N,2:M)+2; %启动值
nn=0;
while max(max(abs(u-uold)))>=10^-10 
        %由于u矩阵为二维，故矩阵最大值的求解利用了两重max算法
        uold=u;  
        u(2:N,2:M)=0.25*(u(3:N+1,2:M)+u(1:N-1,2:M)+u(2:N,3:M+1)+u(2:N,1:M-1));
        nn=nn+1;
end

%figure; surf(X,Y,u')
%figure; contour(X,Y,u',20)
nn
u_error=u'-Z;
figure; meshc(X,Y,u');
colormap hsv
title('差分法得到的温度场的三维图与等高图')
figure; surf(X,Y,u_error)
title('差分法与解析法的误差')
