我正在尝试使用线的方法和 ode15s 命令来解决此处给出的耦合 PDE 系统 - Solve System of PDEs 。引用变量$u1$和$u2$作为$u$和$v$分别，我已将 x 域离散化为-

$$\frac{\partial^2 u}{\partial x^2} = \frac{u_{i+1}-2u_i+u_{i-1}}{(\Delta x)^2} \hspace{10pt} \forall i = 2,3,4...N-1$$ 假设有 N 个网格点。离散边界条件，我得到- $$u_1 = u_2$$ $$v_1 = 0$$ $$v_{N-1} = v_N$$ $$u_N = 1$$ 初始条件是- $$u = 1 \hspace{10pt} \forall \hspace{2pt} i=1,2,3...N$$ $$v = 0\hspace{10pt} \forall \hspace{2pt} i=1,2,3...N$$ 我将我的解决方案向量打包成表格$[u \hspace{4pt}v]^T$这样它看起来像 -$[u_2\hspace{4pt} u_3 \hspace{4pt}.... \hspace{4pt}u_{N-1},\hspace{4pt} v_2 \hspace{4pt} v_3 ....\hspace{4pt} v_{N-1}]$. 所以总的来说，我正在解决$2(N-2)$我可以手动获取变量和边界处的值。这是我到目前为止编写的代码-

clc
clear all
tspan1 = 0:0.005:2;
xmesh1 = 0:0.005:1;

N = size(xmesh1);
N = N(1,2);

N2 = size(tspan1);
N2 = N2(1,2);

y01 = ones(N-2,1);
y02 = zeros(N-2,1);
y0 = [y01; y02];

[t,y] = ode15s(@pde2, tspan1, y0);
u1 = y(:,1:N-2);
u2 = y(:,N-1:2*(N-2));
v1 = u1(:,1);
v2 = ones(N2,1);
w1 = zeros(N2,1);
w2 = u2(:,N-2);
u1 = [v1 u1 v2];
u2 = [w1 u2 w2];

surf(xmesh1, tspan1, u2)
shading interp

function dydt = pde2(t, y)
    xmesh1 = 0:0.005:1;
    N = size(xmesh1);
    N = N(1,2);
    dx = 0.005;
    dydt = zeros(2*(N-2),1);
    dydt(1) = 0.024*(y(2)-y(1))/(dx*dx)-(exp(5.73*y(1))-exp(-11.47*y(N-1)));
    dydt(N-1) = 0.170*(y(N)-2*y(N-1)+0)/(dx*dx)+(exp(5.73*y(1))-exp(-11.47*y(N-1)));
    for i=2:N-2-1
        dydt(i) = 0.024*(y(i+1)-2*y(i)+y(i-1))/(dx*dx)-(exp(5.73*y(i))-exp(-11.47*y(i+N-2)));
        dydt(i+N-2) = 0.170*(y(N-1+i)-2*y(N-2+i)+y(N-3+i))/(dx*dx)+(exp(5.73*y(i))-exp(-11.47*y(i+N-2)));
    end
    dydt(N-2) = 0.024*(1-2*y(N-2)+y(N-3))/(dx*dx)-(exp(5.73*y(N-2))-exp(-11.47*y(2*(N-2)))); 
    dydt(2*(N-2)) = 0.170*(-y(2*(N-2))+y(2*(N-2)-1))/(dx*dx)+(exp(5.73*y(N-2))-exp(-11.47*y(2*(N-2))));
end

这有点低效，但我认为它应该完成工作。问题是我得到了一个与上面链接中给出的非常不同的表面图$u1$和$u2$（要么$u$和$v$）。为了比较，$u2$已附上。

这是原来的情节$u2$. 这就是我得到的-

形状看起来有点相似，但数值相差甚远。我检查了 pde2 函数中的方程，它们对我来说看起来不错。图和值$u1$也有很大不同。有人可以指出其中的错误吗？提前致谢。

编辑：它与方程的刚度有关吗？已经提到解决方案在很小的时候变化很快$t$. 但是我认为使用僵硬的求解器可以规避这个问题。我可以用这种方法很好地解决没有源项的一维热方程，但是这个特殊的问题让我难以理解。