programme Matlab pour la méthode des différences finies pour y « » + e^y = 0

Question

Considérons le problème de la valeur limite suivante:

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

Je suis curieux de savoir comment Matlab résoudra la méthode des différences finies pour ce problème particulier. Je sais que si nous avons une EDO linéaire, par ex. y'' + (e^x)y = 0, avec les mêmes conditions aux limites, alors le programme est assez simple. Disons que nous utilisons la partition de l'intervalle [0,1] dans 20 sous-intervalles égaux, le code suivant fonctionnera:

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

Cependant, il semble que ma question est très différent de cet exemple simple, parce que nous traitent de systèmes d'équations non linéaires. Comment devrions-nous formuler le code dans ce cas?

Answer 1

Vous aurez besoin d'un solveur itératif. Dans le cas le plus simple, résoudre de façon répétée avec

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

Pour une procédure de type Newton, calculer la prochaine approximation u comme ayant une petite différence y afin que e^u=e^y*e^(u-y)=e^y*(1+(u-y)+..) de telle sorte que l'équation linéarisée à résoudre est

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

qui est

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

Si vous souhaitez utiliser solveurs d'équations différentielles Matlab intégré. Vous pouvez utiliser ode45, bvp4c etc. Votre équation peut être ré écrite comme suite d'équations. Laissez y = x1 et ydot = x2, vous obtiendrez x1dot = x2

x2dot = -e^(x1)

Avec vos conditions aux limites, cela peut être résolu en utilisant [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