求解 PDE 方程组 | 您所在的位置:网站首页 › MATLAB怎么求偏微分方程精确解 › 求解 PDE 方程组 |
编写方程代码 在编写方程代码之前,您需要确保它的形式符合 pdepe 求解器的要求: c(x,t,u,∂u∂x)∂u∂t=x-m∂∂x(xmf(x,t,u,∂u∂x))+s(x,t,u,∂u∂x). 在此形式中,PDE 系数是矩阵值,方程变为 [1001]∂∂t[u1u2]=∂∂x[0.024∂u1∂x0.170∂u2∂x]+[-F(u1-u2)F(u1-u2)]. 因此,方程中的系数的值是 m=0 c(x,t,u,∂u∂x)=[11](仅对角线值) f(x,t,u,∂u∂x)=[0.024∂u1∂x0.170∂u2∂x] s(x,t,u,∂u∂x)=[-F(u1-u2)F(u1-u2)] 现在,您可以创建一个函数以编写方程代码。该函数应具有签名 [c,f,s] = pdefun(x,t,u,dudx): x 是独立的空间变量。 t 是独立的时间变量。 u 是关于 x 和 t 微分的因变量。它是二元素向量,其中 u(1) 是 u1(x,t),u(2) 是 u2(x,t)。 dudx 是偏空间导数 ∂u/∂x。它是二元素向量,其中 dudx(1) 是 ∂u1/∂x,dudx(2) 是 ∂u2/∂x。 输出 c、f 和 s 对应于 pdepe 所需的标准 PDE 形式中的系数。 因此,此示例中的方程可由以下函数表示: function [c,f,s] = pdefun(x,t,u,dudx) c = [1; 1]; f = [0.024; 0.17] .* dudx; y = u(1) - u(2); F = exp(5.73*y)-exp(-11.47*y); s = [-F; F]; end(注意:所有函数都作为局部函数包含在示例的末尾。) 编写初始条件代码接下来,编写一个返回初始条件的函数。初始条件应用在第一个时间值处,并为 x 的任何值提供 u(x,t0) 的值。初始条件的数量必须等于方程的数量,因此对于此问题,有两个初始条件。使用函数签名 u0 = pdeic(x) 编写函数。 初始条件为 u1(x,0)=1, u2(x,0)=0. 对应的函数是 function u0 = pdeic(x) u0 = [1; 0]; end 编写边界条件代码现在,编写计算以下边界条件的函数 ∂∂xu1(0,t)=0, u2(0,t)=0,∂∂xu2(1,t)=0, u1(1,t)=1. 对于在区间 a≤x≤b 上提出的问题,边界条件应用于所有 t 以及 x=a 或 x=b。求解器所需的边界条件的标准形式是 p(x,t,u)+q(x,t)f(x,t,u,∂u∂x)=0. 以这种形式编写,u 的偏导数的边界条件需要用通量 f(x,t,u,∂u∂x) 来表示。因此,此问题的边界条件是 对于 x=0,方程为 [0u2]+[10]⋅[0.024∂u1∂x0.170∂u2∂x]=0.
系数是: pL(x,t,u)=[0u2],
qL(x,t)=[10]. 同样,对于 x=1,方程是 [u1-10]+[01]⋅[0.024∂u1∂x0.170∂u2∂x]=0. 系数是: pR(x,t,u)=[u1-10],
qR(x,t)=[01]. 边界函数应使用函数签名 [pl,ql,pr,qr] = pdebc(xl,ul,xr,ur,t): 对于左边界,输入 xl 和 ul 对应于 u 和 x。 对于右边界,输入 xr 和 ur 对应于 u 和 x。 t 是独立的时间变量。 对于左边界,输出 pl 和 ql 对应于 pL(x,t,u) 和 qL(x,t)(对于此问题,x=0)。 对于右边界,输出 pr 和 qr 对应于 pR(x,t,u) 和 qR(x,t)(对于此问题,x=1)。 此示例中的边界条件由以下函数表示: function [pl,ql,pr,qr] = pdebc(xl,ul,xr,ur,t) pl = [0; ul(2)]; ql = [1; 0]; pr = [ur(1)-1; 0]; qr = [0; 1]; end 选择解网格当 t 较小时,此问题的解会快速变化。虽然 pdepe 选择了适合解析急剧变化的时间步,但要在输出绘图中显示该行为,您需要选择适当的输出时间。对于空间网格,在 0≤x≤1 两端的解中都存在边界层,因此您需要在那里指定网格点来解析急剧变化。 x = [0 0.005 0.01 0.05 0.1 0.2 0.5 0.7 0.9 0.95 0.99 0.995 1]; t = [0 0.005 0.01 0.05 0.1 0.5 1 1.5 2];求解方程最后,使用对称性值 m、PDE 方程、初始条件、边界条件以及 x 和 t 的网格来求解方程。 m = 0; sol = pdepe(m,@pdefun,@pdeic,@pdebc,x,t);pdepe 以三维数组 sol 形式返回解,其中 sol(i,j,k) 是在 t(i) 和 x(j) 处计算的解 uk 的第 k 个分量的逼近值。将每个解分量提取到一个单独变量中。 u1 = sol(:,:,1); u2 = sol(:,:,2);对解进行绘图创建在 x 和 t 的所选网格点上绘制的 u1 和 u2 的解的曲面图。 surf(x,t,u1) title('u_1(x,t)') xlabel('Distance x') ylabel('Time t')surf(x,t,u2) title('u_2(x,t)') xlabel('Distance x') ylabel('Time t')局部函数此处列出 PDE 求解器 pdepe 为计算解而调用的局部辅助函数。您也可以将这些函数作为它们自己的文件保存在 MATLAB 路径上的目录中。 function [c,f,s] = pdefun(x,t,u,dudx) % Equation to solve c = [1; 1]; f = [0.024; 0.17] .* dudx; y = u(1) - u(2); F = exp(5.73*y)-exp(-11.47*y); s = [-F; F]; end % --------------------------------------------- function u0 = pdeic(x) % Initial Conditions u0 = [1; 0]; end % --------------------------------------------- function [pl,ql,pr,qr] = pdebc(xl,ul,xr,ur,t) % Boundary Conditions pl = [0; ul(2)]; ql = [1; 0]; pr = [ur(1)-1; 0]; qr = [0; 1]; end % --------------------------------------------- |
CopyRight 2018-2019 实验室设备网 版权所有 |