目的：尝试使用最简单的方法去求解这三个最基本的偏微分方程，为方便绘图，使用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.......

终于发完了、、、

楼主辛苦了，虽然看不懂，但是多谢分享。