ODE BVP: y''-exp(y)=-1
in [Matlab]
|
From: Greg Heath on 30 Dec 2009 01:29 On Dec 24, 3:49 pm, Mirko <mirko.vuko...(a)gmail.com> wrote: > I am trying to solve the following problem: > > y''-exp(y)=-1 > > y(x->-infty)=-1-x^2/2 > y(x->+infty)=0 > > The domain of y is negative. > > I have solved it already one way, discussed below, but > I am wondering if I can solve it as a BVP. > > The way I solved it was to multiply the equation by y', > integrate and obtain > y'=-sqrt(2*(exp(y)-y-1)) > Then rewrite as > dx/sqrt(...) = dy > and integrate. > > Out of curiosity, I would like to know if it is possible > to solve the BVP explicitly. Since the problem is > nonlinear, I have re-written it as > > y''-[exp(y)/y] * y =-1 > > discretized the second order derivative using center > differences, and tried to solve it iteratively. > > But, my first iteration produces values of y>0, which > on the second iteration lead to oscillatory solutions. > > I tried to remove this behavior by using a rescaling factor > y_new = (1-alpha)y_old + alpha* y_new, > > but without success. I am getting some flaky behavior > at the upper boundary. y''- exp(y) = -1 y( x --> -inf ) --> -1 - x^2/2 y( x --> +inf ) --> 0 1. y = 0 solves the ODE and the RHBC 2. From the LHBC y'( x --> -inf) --> -x y''( x --> -inf) --> -1 However, from the ODE, y''( x --> -inf) --> -1 + exp(-1-x^2/2) Therefore, define the asymptotic region x <= x0 < 0 via exp(-1-x0^2/2) = eps where eps is the MATLAB machine epsilon . Therefore eps = 2^(-52) % eps = 2.220446049250313e-16 x0 = -sqrt(-2-2*log(eps)) % x0 = - 8.37181621741867 v0 = y'(x0) = -x0 % v0 = 8.37181621741867 y0 = -1-x0^2/2 % y0 = -36.04365338911716 %%%%%%%%%%%%%%% SYMBOLIC SOLUTION %%%%%%%%%%%%%%%%%%%%%% Integrating the ODE from the LHBC yields (y')^2 = v0^2 + 2*(exp(y)-exp(y0))-2*(y-y0) Integrating again yields the two implicit solutions x-x0 = +/- Int(1/(v0^2+2*(exp(t)-exp(y0))-2*(t-y0)^(1/2),t,y0,y) Only the plus sign yields y' >= 0 when x --> -inf. A complete solution is obtained when the point x = x1 > x0 is found where the explicit solution y = 0 for x >= x1 connects with the implicit solution for x <= x1. IT'S NOT CLEAR HOW TO DO THIS! Note that since the solutions are implicit, BCs cannot be used explicitly in the call of MAPLE's dsolve. For example, >> y = dsolve('D2y = exp(y)-1','x') Warning: Explicit solution could not be found; implicit solution returned. > In C:\MATLAB6p5p1\toolbox\symbolic\dsolve.m at line 292 y = [ Int(1/(2*exp(a)-2*a+C1)^(1/2),a = .. y)-x-C2 = 0] [ Int(-1/(2*exp(a)-2*a+C1)^(1/2),a = .. y)-x-C2 = 0] and >> y=dsolve('(Dy)^2=v0^2+2*(exp(y)-exp(y0))-2*(y-y0)','x') Warning: Explicit solution could not be found; implicit solution returned. > In C:\MATLAB6p5p1\toolbox\symbolic\dsolve.m at line 292 y = [y = log(-lambertw(-exp(1/2*v0^2-exp(y0)+y0)))] [x-Int(1/(v0^2+2*exp(a)-2*exp(y0)-2*a+2*y0)^(1/2),a=.. y)-C1 = 0] [x-Int(-1/(v0^2+2*exp(a)-2*exp(y0)-2*a+2*y0)^(1/2),a=.. y)-C1 = 0] Notice that the solution containing lambertw is independent of x but is nonzero. WHAT IN THE WORLD DOES THAT MEAN? Hope this helps. Greg
From: Greg Heath on 30 Dec 2009 02:16 On Dec 30, 1:29 am, Greg Heath <he...(a)alumni.brown.edu> wrote: > On Dec 24, 3:49 pm, Mirko <mirko.vuko...(a)gmail.com> wrote: > > I am trying to solve the following problem: > > > y''-exp(y)=-1 > > > y(x->-infty)=-1-x^2/2 > > y(x->+infty)=0 > > > The domain of y is negative. > > > I have solved it already one way, discussed below, but > > I am wondering if I can solve it as a BVP. > > > The way I solved it was to multiply the equation by y', > > integrate and obtain > > y'=-sqrt(2*(exp(y)-y-1)) > > Then rewrite as > > dx/sqrt(...) = dy > > and integrate. > > > Out of curiosity, I would like to know if it is possible > > to solve the BVP explicitly. Since the problem is > > nonlinear, I have re-written it as > > > y''-[exp(y)/y] * y =-1 > > > discretized the second order derivative using center > > differences, and tried to solve it iteratively. > > > But, my first iteration produces values of y>0, which > > on the second iteration lead to oscillatory solutions. > > > I tried to remove this behavior by using a rescaling factor > > y_new = (1-alpha)y_old + alpha* y_new, > > > but without success. I am getting some flaky behavior > > at the upper boundary. > > y''- exp(y) = -1 > > y( x --> -inf ) --> -1 - x^2/2 > y( x --> +inf ) --> 0 > > 1. y = 0 solves the ODE and the RHBC > 2. From the LHBC > > y'( x --> -inf) --> -x > y''( x --> -inf) --> -1 > > However, from the ODE, > > y''( x --> -inf) --> -1 + exp(-1-x^2/2) > > Therefore, define the asymptotic region x <= x0 < 0 via > > exp(-1-x0^2/2) = eps > > where eps is the MATLAB machine epsilon . Therefore > > eps = 2^(-52) % eps = 2.220446049250313e-16 > x0 = -sqrt(-2-2*log(eps)) % x0 = - 8.37181621741867 > v0 = y'(x0) = -x0 % v0 = 8.37181621741867 > y0 = -1-x0^2/2 % y0 = -36.04365338911716 > > %%%%%%%%%%%%%%% SYMBOLIC SOLUTION %%%%%%%%%%%%%%%%%%%%% -----SNIP %%%%%%%%%%%%%%% NUMERICAL SOLUTION %%%%%%%%%%%%%%%%%%%%%% close all, clear all, clc, k = 0; x0 = -sqrt(-2*log(eps)-2) % -8.3718 dx = 0.05 N = 491 % For a symmetrical plot x = zeros(N,1); v = zeros(N,1); % v = y' y = zeros(N,1); x(1) = x0-dx; % -8.4218 v(1) = -(x0-dx); % 8.4218 y(1) = -1-(x0-dx)^2/2; % -36.463 x(2) = x0; v(2) = -x0; y(2) = -1-x0^2/2; % -36.044 for i = 2:N-1 x(i+1) = x(i)+dx; v(i+1) = v(i-1) + 2*dx*(exp(y(i))-1); y(i+1) = y(i-1) + 2*dx*v(i); end k=k+1; figure(k) plot(x,y,'.') [ymax imax] = max(y); % imax = 246 itab = [1 imax N]'; summary = [ itab x(itab) y(itab)] % SUMMARY % % index x y % 1 -8.4218 -36.463 % 246 3.8282 -0.028995 % ymax % 491 16.078 -36.465 % COMMENTS % % The plot looks symmetric about x1 = x(246) % with a shape something like y = ymax-a*(x-x1)^2 % except that it is much flatter near ymax. % % I suspect that if I had chosen x0 more negative, % ymax = 0 and that would be the point where this % numerical solution connects with the symbolic % solution y = 0 for x >= x1. % % Hope this helps. % % Greg
From: Greg Heath on 30 Dec 2009 04:46 On Dec 30, 2:16 am, Greg Heath <he...(a)alumni.brown.edu> wrote: > On Dec 30, 1:29 am, Greg Heath <he...(a)alumni.brown.edu> wrote: > > On Dec 24, 3:49 pm, Mirko <mirko.vuko...(a)gmail.com> wrote: > > > I am trying to solve the following problem: > > > > y''-exp(y)=-1 > > > > y(x->-infty)=-1-x^2/2 > > > y(x->+infty)=0 > > > > The domain of y is negative. > > > > I have solved it already one way, discussed below, but > > > I am wondering if I can solve it as a BVP. > > > > The way I solved it was to multiply the equation by y', > > > integrate and obtain > > > y'=-sqrt(2*(exp(y)-y-1)) > > > Then rewrite as > > > dx/sqrt(...) = dy > > > and integrate. > > > > Out of curiosity, I would like to know if it is possible > > > to solve the BVP explicitly. Since the problem is > > > nonlinear, I have re-written it as > > > > y''-[exp(y)/y] * y =-1 > > > > discretized the second order derivative using center > > > differences, and tried to solve it iteratively. > > > > But, my first iteration produces values of y>0, which > > > on the second iteration lead to oscillatory solutions. > > > > I tried to remove this behavior by using a rescaling factor > > > y_new = (1-alpha)y_old + alpha* y_new, > > > > but without success. I am getting some flaky behavior > > > at the upper boundary. > > > y''- exp(y) = -1 > > > y( x --> -inf ) --> -1 - x^2/2 > > y( x --> +inf ) --> 0 > > > 1. y = 0 solves the ODE and the RHBC > > 2. From the LHBC > > > y'( x --> -inf) --> -x > > y''( x --> -inf) --> -1 > > > However, from the ODE, > > > y''( x --> -inf) --> -1 + exp(-1-x^2/2) > > > Therefore, define the asymptotic region x <= x0 < 0 via > > > exp(-1-x0^2/2) = eps > > > where eps is the MATLAB machine epsilon . Therefore > > > eps = 2^(-52) % eps = 2.220446049250313e-16 > > x0 = -sqrt(-2-2*log(eps)) % x0 = - 8.37181621741867 > > v0 = y'(x0) = -x0 % v0 = 8.37181621741867 > > y0 = -1-x0^2/2 % y0 = -36.04365338911716 > > > %%%%%%%%%%%%%%% SYMBOLIC SOLUTION %%%%%%%%%%%%%%%%%%%%% > > -----SNIP > %%%%%%%%%%%%%%% NUMERICAL SOLUTION %%%%%%%%%%%%%%%%%%%%%% > > close all, clear all, clc, k = 0; > > x0 = -sqrt(-2*log(eps)-2) % -8.3718 > dx = 0.05 > N = 491 % For a symmetrical plot > x = zeros(N,1); > v = zeros(N,1); % v = y' > y = zeros(N,1); > > x(1) = x0-dx; % -8.4218 > v(1) = -(x0-dx); % 8.4218 > y(1) = -1-(x0-dx)^2/2; % -36.463 > x(2) = x0; > v(2) = -x0; > y(2) = -1-x0^2/2; % -36.044 > > for i = 2:N-1 > x(i+1) = x(i)+dx; > v(i+1) = v(i-1) + 2*dx*(exp(y(i))-1); > y(i+1) = y(i-1) + 2*dx*v(i); > end > > k=k+1; figure(k) > plot(x,y,'.') > > [ymax imax] = max(y); % imax = 246 -----SNIP [absvmin jmin] = min(abs(v)); % jmin = 246 itab = [1 imax N]'; summary = [ itab x(itab) v(itab) y(itab)] % SUMMARY % % index x v y % 1 -8.422 8.422 -36.463 % 246 3.828 -2.163e-6 -0.029 % 491 16.078 -8.422 -36.465 % COMMENTS % % The plot looks symmetric about x1 = 3.83 with % % min(abs(y')) = abs(y'(x1)) = 2.16e-6 % ymax = y(x1) = -0.029 % % and a shape something like y = ymax-a*(x-x1)^2 % except that it is much flatter near ymax. % % I suspect that if I had chosen x0 more negative, % ymax = 0 and that would be the point where this % numerical solution connects with the symbolic % solution y = 0 for x >= x1. Maybe not! Replacing eps with eps^2 and using N = 633 yields index x v y 1 -11.974 11.974 -72.685 317 3.826 5.164e-5 -0.029 Same ymax! 633 19.626 -11.97 -72.641 Although x1 = 3.83, ymax =y(x1) = -0.029 has not changed abs(y'(x1)) increases to 5.16e-5. A more accurate finite difference solution would yield more confidence to my conjecture. Hope this helps. Greg
From: Greg Heath on 30 Dec 2009 18:02 On Dec 24, 3:49 pm, Mirko <mirko.vuko...(a)gmail.com> wrote: > Hello, > > I am trying to solve the following problem: > > y''-exp(y)=-1 > > y(x->-infty)=-1-x^2/2 > y(x->+infty)=0 > > The domain of y is negative. > > I have solved it already one way, discussed below, but I am wondering > if I can solve it as a BVP. > > The way I solved it was to multiply the equation by y', integrate and > obtain > y'=-sqrt(2*(exp(y)-y-1)) > Then rewrite as > dx/sqrt(...) = dy > and integrate. > > Out of curiosity, I would like to know if it is possible to solve the > BVP explicitly. Since the problem is nonlinear, I have re-written it > as > > y''-[exp(y)/y] * y =-1 > > discretized the second order derivative using center differences, and > tried to solve it iteratively. > > But, my first iteration produces values of y>0, which on the second > iteration lead to oscillatory solutions. > > I tried to remove this behavior by using a rescaling factor > y_new = (1-alpha)y_old + alpha* y_new, > > but without success. I am getting some flaky behavior at the upper > boundary. There are 3 solution branches: y > 0, y= 0 and y < 0. If you start iterations near y = 0 there is a risk of jumping from one branch to another, causing nonconvergent oscillations. This might be able to be mitigated by using a smaller deltax. Howevr, it is obviously better to start iterations in the x --> -inf asymptotic region y ~ -1 - x^2/2 < 0, y' ~ -x > 0 y'' ~ -1 < 0 and match up with y = 0 solution to obey the x --> inf BC. Hope this helps. Greg.
From: Greg Heath on 1 Jan 2010 14:25
On Dec 30 2009, 4:46 am, Greg Heath <he...(a)alumni.brown.edu> wrote: > On Dec 30, 2:16 am, Greg Heath <he...(a)alumni.brown.edu> wrote: > > On Dec 30, 1:29 am, Greg Heath <he...(a)alumni.brown.edu> wrote: > > > On Dec 24, 3:49 pm, Mirko <mirko.vuko...(a)gmail.com> wrote: > > > > I am trying to solve the following problem: > > > > > y''-exp(y)=-1 > > > > > y(x->-infty)=-1-x^2/2 > > > > y(x->+infty)=0 -----SNIP > % COMMENTS > % > % The plot looks symmetric about x1 = 3.83 with > % > % min(abs(y')) = abs(y'(x1)) = 2.16e-6 > % ymax = y(x1) = -0.029 > % > % and a shape something like y = ymax-a*(x-x1)^2 > % except that it is much flatter near ymax. > % > % I suspect that if I had chosen x0 more negative, > % ymax = 0 and that would be the point where this > % numerical solution connects with the symbolic > % solution y = 0 for x >= x1. > > Maybe not! > > Replacing eps with eps^2 and using N = 633 yields > > index x v y > 1 -11.974 11.974 -72.685 > 317 3.826 5.164e-5 -0.029 Same ymax! > 633 19.626 -11.97 -72.641 > > Although x1 = 3.83, ymax =y(x1) = -0.029 > has not changed abs(y'(x1)) increases to > 5.16e-5. > > A more accurate finite difference solution would > yield more confidence to my conjecture. Or significantly refute it! Instead of decreasing negative x0, I (duh) decreased dx. close all, clear all, clc, k = 0; format long g % Previous runs % % dx = 0.05 % x0 = -sqrt(-2*log(eps)-2) % -8.372 % N = 491 % Chosen for plot symmetry % % x0 = -sqrt(-2*log(eps^2)-2) % -11.974 % N = 633 % % x0 = -sqrt(-2*log(eps^3)-2) % -14.688 % N = 741 dx = 0.0001 x0 = -2 N = 68687 x = zeros(N,1); v = zeros(N,1); % v = y' y = zeros(N,1); x(1) = x0-dx; % v(1) = -(x0-dx); % y(1) = -1-(x0-dx)^2/2; % x(2) = x0; v(2) = -x0; y(2) = -1-x0^2/2; % for i = 2:N-1 x(i+1) = x(i)+dx; v(i+1) = v(i-1) + 2*dx*(exp(y(i))-1); y(i+1) = y(i-1) + 2*dx*v(i); end [ymax imax] = max(y) % [ -0.333, 34344] x1 = x(imax) % 1.434 [ymin imin] = min(y) % [ -3.001, 68687 ] x2 = x(imin) % 4.8769 [absvmin jmin]=min(abs(v)) % [2.32e-5,34345] itab = [1 imax jmin N]'; summary = [ itab x(itab) v(itab) y(itab)] % index x v y % 1 -2.000 2.000 -3.000 % 34344 1.434 -2.354e-005 -0.333 % 34345 1.434 -2.316e-005 -0.333 % 68687 4.869 -2.000 -3.000 yinf = -1-x.^2/2; xinf = x(yinf >= ymin); yinf = yinf(yinf >= ymin); k=k+1; figure(k), hold on % figure 1 plot(x,zeros(N,1),'k') plot(x,y) plot(xinf,yinf,'r') xlim([x0-1 x(end)+1]) ylim([ymin-1 ymax+1]) % COMMENTS % % 1. The y plot looks like it is symmetric about % x = x1 = 1.434 where ymax = -0.333 and y' ~ 2e-5. % 2. A stepwise polynomial fit to the numerical result % yields % % y = -0.328 - 0.150*(x-x1)^2 - 0.00661*(x-x1)^4 % % with R^2 = 0.9996. % 3. Most importantly, ymax is so far away from 0 % that I have to retract my conjecture that with % increased precision ymax = 0 and the negative % numerical solution can be joined to the y = 0 % theoretical result that satisfies the RHBC. BOTTOM LINE: I cannot find a negative y solution that will connect to the y = 0 solution that satisfies the RHBC. HELP! Greg |