программа MATLAB для метода конечных разностей для у «» + е^у = 0

Question

Рассмотрим следующую краевую задачу:

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

Мне любопытно, как MATLAB будет решить метод конечных разностей для этой конкретной проблемы. Я знаю, что если у нас есть линейное ОДУ, т.е. 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'; 

Однако, кажется, что мой вопрос очень отличается от этого простого примера, потому что мы имеют дело с системами нелинейных уравнений. Как мы должны формулировать код в этом случае?

Answer 1

Вам понадобится итеративный решатель. В самом простом случае многократно решить с

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

Answer 2

Если вы хотите использовать 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