有限差分法求解一维导热方程、一维波动方程、二维拉普拉斯方程

目的：尝试使用最简单的方法去求解这三个最基本的偏微分方程，为方便绘图，使用matlab软件。
求解过程：1、首先写出方程和定解条件
                  2、其次将方程最简单的离散化，并改写为矢量形式，并给定边界条件和初值
                  3、再编程求解，绘图得到结果

一维非稳态导热方程：1-3部分在图片1。4.matlab代码如下：
%% 给定边值，时间间隔，空间间隔,稳定性条件
u_a=1;    %% 左边值
u_b=1;    %% 右边值
deta_t=0.0001;  %% 时间间隔
Total_time=0.5; %% 总时长
deta_x=0.02;    %% 空间间隔
n=1:51;         %% 节点序号
sigma=deta_t/(deta_x^2);
k=0;            %% 时刻标志位，用于绘图
%% 一维杆的初值
u(1)=u_a;       %%  左边值
u(51)=u_b;      %%  右边值
for i=2:50
    u(i)=(i-1)/50;    %% 初始时刻的函数
end
plot((n-1)/50,u,'r');   %% 画出初始时刻温度图
hold on
%% 时间步进
for j=0:deta_t:Total_time
    u_next(1)=u(1);     %% 左边界值，在该问题中不随时间变化
    u_next(51)=u(51);   %% 左边界值，在该问题中不随时间变化
    %% 计算下一时刻温度值
    for i=2:50
        u_next(i)=sigma*u(i-1)+(1-2*sigma)*u(i)+sigma*u(i+1);
    end
    %%  交换数据与绘图
    u=u_next;
    k=k+1;
    if(mod(k,1000)==0)
        plot((n-1)/50,u,'-');
    end
end
%% 绘图修正
xlabel('X值');
ylabel('温度');
legend('0秒','0.1秒','0.2秒','0.3秒','0.4秒','0.5秒');
axis([0 1 0 1.2]);



5. 结果分析：随着时间的增长，左右边界为温度逐渐向内部传达，最终温度值将达到统一，符合热力学第二定律。

一维波动方程：1-3部分在图片1。
4.matlab代码如下：

%% 给定边值，时间间隔，空间间隔,稳定性条件u_a=1;    %% 左边值
u_b=1;    %% 右边值deta_t=0.001;  %% 时间间隔
Total_time=0.25; %% 总时长deta_x=0.02;    %% 空间间隔
n=1:51;         %% 节点序号nn=max(n-1);    %% 最大步数减一
sigma=(deta_t^2)/(deta_x^2);k=0;            %% 时刻标志位，用于绘图
%% 一维杆的初值u(1)=u_a;        %%  左边值
u(nn+1)=u_b;     %%  右边值for i=2:nn
   u(i)=cos(2*pi*(i-1)/nn);    %% 初始时刻的函数end
plot((n-1)/nn,u,'ro');   %% 画出初始时刻位移图hold on
%%  起始步u_next(1)=u(1);     %% 左边界值，在该问题中不随时间变化
u_next(nn+1)=u(nn+1);   %% 左边界值，在该问题中不随时间变化    %% 计算第一时刻位置值
   for i=2:nn        u_next(i)=0.5*(sigma*u(i+1)+(2-2*sigma)*u(i)+sigma*u(i-1));
   end    u_last=u;
   u=u_next;    k=k+1;
plot((n-1)/nn,u,'b+');   %% 画出起始步时刻位移图%% 时间步进
for j=deta_t:deta_t:Total_time    u_next(1)=u(1);     %% 左边界值，在该问题中不随时间变化
   u_next(nn+1)=u(nn+1);   %% 左边界值，在该问题中不随时间变化    %% 计算下一时刻位移值
   for i=2:nn        u_next(i)=sigma*u(i+1)+(2-2*sigma)*u(i)+sigma*u(i-1)-u_last(i);
   end    %%  交换数据与绘图
   u_last=u;    u=u_next;
   k=k+1;    if(mod(k,50)==0)
       plot((n-1)/nn,u,'-');    end
end%% 绘图修正
xlabel('X值');ylabel('位移');
legend('0秒','0.001秒','0.05秒','0.10秒','0.15秒','0.20秒','0.25秒');axis([0 1 -1 1]);5 结果分析：由于太困了，仅根据改变不同边界条件和初值，曲线表现得像两端捏死的绳子的振动过程，大致认为该程序正确。

如有错误 请告知

二维拉普拉斯方程：1-3部分在图片1、2。
4.matlab代码如下：
%% 给定空间区域，空间间隔,超松弛因子
deta_x=0.01;  %% x间隔
Total_x=1; %% 总x长
deta_y=0.01;    %% y间隔
Total_y=1; %% 总y长
n=1:((Total_x/deta_x)+1);         %% 节点序号
nn=max(n)-1;
sigma=0.25;
omega=2/(1+sqrt(1-
((cos(pi/nn)+cos(pi/nn))/2 ).^2) );  %% 最佳超松弛因子
iteration_max=10000;    %迭代总次数
convergence=1e-8;       %迭代精度
%%  二维温度场的坐标值和初值
for i=1:nn
   x(i)=i*deta_x; %% 定义网格点坐标
end
for j=1:nn
   y(j)=j*deta_y; %% 定义网格点坐标
end
u=zeros(nn,nn);
%对u赋初值
%%  二维温度场的四条边的边界值
for i=1:nn
   j=1;
   u(i,j)=1;       %% y=0的u值
end
for i=1:nn
   j=nn;
   u(i,j)=0;       %% y=1的u值
end
for j=1:nn
   i=1;
   u(i,j)=1;       %% x=0的u值
end
for j=1:nn
   i=nn;
   u(i,j)=1;       %% x=1的u值
end
u_last=u;
u_next=u;
%% 设置迭代
for
iteration=1:iteration_max
   for j=2:nn-1
       for i=2:nn-1
         
u_next(i,j)=omega.*sigma.*(u(i+1,j)+u(i,j+1)+u(i-1,j)+u(i,j-1))...
               +(1-omega)*u(i,j);
           u(i,j)=u_next(i,j);
       end
   end
   if( max(max(u-u_last))< convergence)
       break;
   end
   u_last=u;
end
%% 绘图修正
[C,h]=contourf(x,y,u);
V=0:0.1:1;
clabel(C,h,V,'FontSize',12,'Color','red','LabelSpacing',1200);
xlabel('X值');
ylabel('Y值');
title('等值线云图')
axis([0 1 0 1]);
5.结果分析
容许误差为1e-8,最大迭代次数10000次，实际迭代次数335次。雅各比迭代次数或高斯赛德尔迭代相比较收敛较快。将结果绘制成云图如图3，较符合物理规律。
6.注1：椭圆型方程也常用于有限差分法的网格变换。
注2：初值条件打错了，是x=1时，u=0.......

终于发完了、、、