%example-4,用松弛法（Gauss-Seidel迭代法）求解温度场
clear all
close all
clc

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(1),surf(X,Y,Z)
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); %边界条件
u(2:N,2:M)=u(2:N,2:M)+2.0;                    %初始猜测值
%unew=u;
omega=0.9; %松弛因子，omega=1时为gauss-seidel迭代法  
n=0; %迭代次数
err=1;
    while abs(err)>=10^-10 
        n=n+1;
        err=0;
        for i=2:N
            for j=2:M
                %Gauss-Seidel迭代法
                relx=omega/4*(u(i,j+1)+u(i,j-1)+u(i+1,j)+u(i-1,j)-4*u(i,j));
                u(i,j)=u(i,j)+relx;
                if abs(relx)>err
                    err=abs(relx);
                end
                %err=abs(relx/omega);
            end
        end
    end

n
abs(err)
figure(2); surf(X,Y,u')
title('Gauss-Seidel迭代法得到的温度场的三维图')
u_error=u'-Z;
figure(3); surf(X,Y,u_error)
title('数值方法与解析方法的误差')
