目的:尝试使用最简单的方法去求解这三个最基本的偏微分方程,为方便绘图,使用matlab软件。
求解过程:1、首先写出方程和定解条件 2、其次将方程最简单的离散化,并改写为矢量形式,并给定边界条件和初值 3、再编程求解,绘图得到结果 |
|
|
飞行棋狂魔
发布于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. 结果分析:随着时间的增长,左右边界为温度逐渐向内部传达,最终温度值将达到统一,符合热力学第二定律。 |
|
|
飞行棋狂魔
发布于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 结果分析:由于太困了,仅根据改变不同边界条件和初值,曲线表现得像两端捏死的绳子的振动过程,大致认为该程序正确。 如有错误 请告知
|
|
|
叽叽歪
发布于2018-04-08 12:53
地板F
|
|
|
飞行棋狂魔
发布于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=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....... |
|
|
飞行棋狂魔
发布于2018-04-08 15:35
5楼F
终于发完了、、、
|
|
|