主要内容

方域的波动方程

这个例子演示了如何用solvepde函数。

标准的二阶波动方程是

2 u t 2 - u 0

要用工具箱的形式表达这一点,请注意solvepde函数解决了表单的问题

2 u t 2 - c u + 一个 u f

标准波动方程有系数 1 c 1 一个 0 , f 0

c = 1;= 0;f = 0;m = 1;

在方形域上解决问题。的squareg函数描述这个几何图形。创建一个模型对象,并包括几何体。绘制几何图形并查看边缘标签。

numberOfPDE = 1;模型= createpde (numberOfPDE);geometryFromEdges(模型、@squareg);pdegplot(模型,“EdgeLabels”“上”);ylim ([-1.1 - 1.1]);轴平等的标题“显示边缘标签的几何图形”;包含xylabely

图中包含一个坐标轴。标题为“显示几何边缘标签”的轴包含5个类型为行、文本的对象。

指定PDE系数。

specifyCoefficients(模型,“米”米,' d ',0,“c”c“一个”一个,“f”f);

左边(边4)和右边(边2)的Dirichlet边界条件设为零,上边(边1)和下边(边3)的Neumann边界条件设为零。

applyBoundaryCondition(模型,“边界条件”“边缘”(2、4),“u”, 0);applyBoundaryCondition(模型,“纽曼”“边缘”3 ([1]),‘g’, 0);

创建并查看问题的有限元网格。

generateMesh(模型);图pdemesh(模型);ylim ([-1.1 - 1.1]);轴平等的包含xylabely

图中包含一个坐标轴。轴线包含2个线型对象。

设置初始条件:

  • u x 0 反正切 因为 π x 2

  • u t | t 0 3. π x 经验值 π y 2

U0 = @(location) atan(cos(pi/2*location.x)); / /位置Ut0 = @(location) 3*sin(pi*location.x).*exp(sin(pi/2*location.y));setInitialConditions(模型、情况ut0);

这种选择避免了将能量放入更高的振动模式,并允许一个合理的时间步长。

指定解乘以31个从0到5等间距的时间点。

n = 31;tlist = linspace (0 5 n);

设置SolverOptions。ReportStatistics模型“上”

model.SolverOptions.ReportStatistics =“上”;结果= solvepde(模型、tlist);
459个成功的步骤39个失败的尝试998个函数计算1个偏导数115个LU分解997个线性系统的解金宝搏官方网站
u = result.NodalSolution;

创建一个动画来可视化所有时间步骤的解决方案。保持一个固定的垂直比例尺,首先计算的最大值和最小值u在所有时间,缩放所有的地块来使用这些 z 设在限制。

图umax = max(max(u));umin = min (min (u));I = 1:n pdeplot(模型,“XYData”u(:,我),“ZData”u(:,我),“ZStyle”“连续”“网”“关闭”);轴([-1 1 -1 1 umin umax]);caxis ([umin umax]);包含xylabelyzlabeluM (i) = getframe;结束

图中包含一个坐标轴。坐标轴包含一个patch类型的对象。

要播放动画,请使用电影(M)命令。