• UID60
  • 登录2018-04-08
  • 粉丝26
  • 发帖122
  • 科研点数9点
MT论坛理事会理事
Maker
二级优质勋章
优异服役勋章
论坛认证学者
飞行棋狂魔 发布于2018-02-28 14:01
5/1397

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

楼层直达
目的:尝试使用最简单的方法去求解这三个最基本的偏微分方程,为方便绘图,使用matlab软件。
求解过程:1、首先写出方程和定解条件
                  2、其次将方程最简单的离散化,并改写为矢量形式,并给定边界条件和初值
                  3、再编程求解,绘图得到结果
看完回个贴呗,反正不要钱QAQ
  • UID60
  • 登录2018-04-08
  • 粉丝26
  • 发帖122
  • 科研点数9点
MT论坛理事会理事
Maker
二级优质勋章
优异服役勋章
论坛认证学者
飞行棋狂魔 发布于2018-04-06 18:23
沙发F
一维非稳态导热方程: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. 结果分析:随着时间的增长,左右边界为温度逐渐向内部传达,最终温度值将达到统一,符合热力学第二定律。
看完回个贴呗,反正不要钱QAQ
  • UID60
  • 登录2018-04-08
  • 粉丝26
  • 发帖122
  • 科研点数9点
MT论坛理事会理事
Maker
二级优质勋章
优异服役勋章
论坛认证学者
飞行棋狂魔 发布于2018-04-07 02:33
板凳F
一维波动方程: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 结果分析:由于太困了,仅根据改变不同边界条件和初值,曲线表现得像两端捏死的绳子的振动过程,大致认为该程序正确。
如有错误 请告知

看完回个贴呗,反正不要钱QAQ
  • UID6
  • 登录2018-08-25
  • 粉丝86
  • 发帖274
  • 科研点数4点
高级电子或计算机学者
2016年论坛年会举办者
2018年论坛年会举办者
MT论坛理事会理事
论坛之星服役勋章
追星Maker纪念章
众望瞩目勋章
荣誉会员勋章
优异服役勋章
Maker
叽叽歪 发布于2018-04-08 12:53
地板F
站长反馈邮箱:YBA@makertime.org
  • UID60
  • 登录2018-04-08
  • 粉丝26
  • 发帖122
  • 科研点数9点
MT论坛理事会理事
Maker
二级优质勋章
优异服役勋章
论坛认证学者
飞行棋狂魔 发布于2018-04-08 15:28
4楼F
二维拉普拉斯方程: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=0u
end
for i=1:nn
   j=nn;
   u(i,j)=0;       %% y=1u
end
for j=1:nn
   i=1;
   u(i,j)=1;       %% x=0u
end
for j=1:nn
   i=nn;
   u(i,j)=1;       %% x=1u
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.......
看完回个贴呗,反正不要钱QAQ
  • UID60
  • 登录2018-04-08
  • 粉丝26
  • 发帖122
  • 科研点数9点
MT论坛理事会理事
Maker
二级优质勋章
优异服役勋章
论坛认证学者
飞行棋狂魔 发布于2018-04-08 15:35
5楼F
终于发完了、、、
看完回个贴呗,反正不要钱QAQ
您需要登录后才可以回帖
发表回复

杩斿洖椤堕儴