y'' + e^y = 0 的有限差分法的 MATLAB 程序答案

【问题标题】：MATLAB program for finite difference method for y'' + e^y = 0y'' + e^y = 0 的有限差分法的 MATLAB 程序
【发布时间】：2017-01-12 22:13:21
【问题描述】：

考虑以下边值问题：

y'' + e^y = 0 i.e. y(0) = y(1) = 0.

我很好奇 MATLAB 将如何解决这个特定问题的有限差分法。我知道如果我们有一个线性 ODE，例如y'' + (e^x)y = 0，同样的边界条件，那么程序就相当简单了。假设我们将区间 [0,1] 划分为 20 个相等的子区间，那么下面的代码将起作用：

N = 19;
h = 1/N;
x = linspace(0, 1, N+1)';
A(1,1) = 1;
F(1) = 0;

for k=2:N
  A(k,k-1) = -1/h^2;
  A(k,k)   = 2/h^2+exp(x(k));
  A(k,k+1) = -1/h^2;
  F(k) = 0;
end

A(N+1, N+1)=1;
F(N+1) = 0;
U = A\F';

但是，我的问题似乎与这个简单的例子有很大不同，因为我们正在处理非线性方程组。在这种情况下我们应该如何制定代码？

【问题讨论】：

标签： matlab solver differential-equations

【解决方案1】：

如果您想使用 Matlab 内置的微分方程求解器。您可以使用ode45, bvp4c 等。您的方程可以重写为以下方程组。让y = x1和ydot = x2，你会得到 x1dot = x2

x2dot = -e^(x1)

根据您的边界条件，可以使用[bvp4c]1 解决此问题

function SOQ
solinit = bvpinit(linspace(0,1,5),[0 0]);% initial guess taken as [0 0]
sol = bvp4c(@ode,@bouncond,solinit);
x = linspace(0,1);
y = deval(sol,x);
plot(x,y(1,:));
end

function dydx = ode(x,y) % system of equations
dydx = [y(2);-exp(y(1))];
end

function res = bouncond(ya,yb) % boundary conditions
res = [ya(1);yb(1)];
end

【讨论】：

【解决方案2】：

您将需要一个迭代求解器。在最简单的情况下反复求解

  A(k,k-1) = -1/h^2;
  A(k,k)   = 2/h^2;
  A(k,k+1) = -1/h^2;

  F(k) = -exp(y(k));

对于类似牛顿的过程，计算下一个近似值 u 与 y 的差异很小，因此 e^u=e^y*e^(u-y)=e^y*(1+(u-y)+..) 使得要求解的线性化方程为

u'' + e^y*u = F(x) = -e^y*(1-y)

那是

  A(k,k-1) = -1/h^2;
  A(k,k)   = 2/h^2 + exp(y(k));
  A(k,k+1) = -1/h^2;

  F(k) = -exp(y(k))*(1-y(k));

【讨论】：