有限差分法求解一维导热方程、一维波动方程、二维拉普拉
目的：尝试使用最简单的方法去求解这三个最基本的偏微分方程，为方便绘图，使用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();
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();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
%% 绘图修正
=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();
5.结果分析
容许误差为1e-8,最大迭代次数10000次，实际迭代次数335次。雅各比迭代次数或高斯赛德尔迭代相比较收敛较快。将结果绘制成云图如图3，较符合物理规律。
6.注1：椭圆型方程也常用于有限差分法的网格变换。
注2：初值条件打错了，是x=1时，u=0....... 终于发完了、、、 楼主辛苦了，虽然看不懂，但是多谢分享。:victory:
页:
[1]