主要内容

方域波动方程

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

标准的二阶波动方程是

2 u t 2 - u 0

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

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

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

C = 1;A = 0;F = 0;M = 1;

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

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

图中包含一个轴对象。标题为Geometry with Edge Labels的axis对象包含5个类型为line, text的对象。

指定PDE系数。

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

在左侧(边4)和右侧(边2)设置零狄利克雷边界条件,在顶部(边1)和底部(边3)设置零诺伊曼边界条件。

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

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

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

图中包含一个轴对象。axis对象包含2个line类型的对象。

设置以下初始条件:

  • 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);

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

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

N = 31;Tlist = linspace(0,5,n);

设置SolverOptions。ReportStatistics模型“上”

model.SolverOptions.ReportStatistics =“上”;Result = solvepde(model,tlist);
457个成功步骤38个失败尝试992个函数求值1个偏导数113个LU分解991个线性系统解
u = result. nodesolution;

创建一个动画来可视化所有时间步骤的解决方案。的最大值和最小值,保持一个固定的垂直刻度u在所有的时间里,缩放所有的地块来使用它们 z 设在限制。

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

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

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

Baidu
map