Prev: How to get Fourier descriptor
Next: (Bode)2 figure in one
From: Greg Heath on 1 Jan 2010 17:18
On Dec 30 2009, 1:29 am, Greg Heath <he...(a)alumni.brown.edu> wrote:
>-----SNIP
> 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!

Indeed. I have previous posts that not only failed,
but they left me with no other ideas to try.

> Note that since the solutions are implicit, BCs
> cannot be used explicitly in the call of MAPLE's
> dsolve.
>
> For example,

-----SNIP
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?

The solution containing lambertw is the explicit solution for
(Dy)^2 == 0:

0 = v0^2+2*(exp(y)-exp(y0))-2*(y-y0)

y - exp(y) = v0^2/2 + y0 - exp(y0)
exp(y - exp(y)) = exp(v0^2/2 + y0-exp(y0))
-exp(y)*exp(-exp(y)) = -exp(v0^2/2+y0 -exp(y0))

which has the form of lambert's equation solved
by

exp(y) = -lambertw(-exp(v0^2/2+y0 -exp(y0)))

which, in turn can be solved by taking the log.

However, note that direct substitution into MAPLE
yields error messages for:

y = dsolve('0=v0^2+2*(exp(y)-exp(y0))-2*(y-y0)','x')

(It's no longer an ODE) and

y = solve('0=v0^2+2*(exp(y)-exp(y0))-2*(y-y0)','x')

(The equation doesn't contain x). Whereas

clear x y v0 y0
solve('0=v0^2+2*(exp(y)-exp(y0))-2*(y-y0)')

surprisingly(No declared dependent variable),
yields

y = -lambertw(-exp(1/2*v0^2-exp(y0)+y0))
+1/2*v0^2-exp(y0)+y0

which can be shown to solve the original equation.

BOTTOM LINE:

(Dy)^2 == 0 ==> v0 = 0

==> y - exp(y) = y0 - exp(y0) ==> y = y0

Hope this helps.

Greg

From: Greg Heath on 3 Jan 2010 10:13
On Dec 24, 3:49 pm, Mirko <mirko.vuko...(a)gmail.com> wrote:

-----SNIP

Numerical Solution using MATLAB

% y'' - exp(y) = -1 % ODE
% y( x--> -inf) --> -1- x^2/2 % LHBC
% y( x--> +inf) --> 0 % RHBC
% From LHBC:
% y'( x--> -inf) --> - x
% y''( x--> -inf) --> - 1
% Combining ODE and LHBC
% y''(x--> -inf) --> -1 + exp(-1- x^2/2)
% Therefore define LHBC asymptotic region by
% exp(-1-x0^2/2) << 1, x0 < 0
% For example,
% x0 <= -2 ==> exp(-1-x0^2/2) < 0.05
% From ODE and RHBC
% y''(x--> +inf) --> -1 + (1 + y)
% Therefore
% y(x --> +inf) --> A*exp(-x+x2), unknown A
% y'(x --> + inf) --> -y
% Define
% v = y'
% v(x --> -inf) --> -x
% v(x --> +inf) --> -y

close all, clear all, clc, k = 0;
format short g

% THE LEFT HAND BOUNDARY CONDITION SOLUTION

x0 = -2
dx = 0.01
N = 687 % Chosen for plot symmetry
x = x0+dx*(0:N-1)';

% LHBC
yLinf = -1-x.^2/2; % LH asymptote

%Initialization
vL = -x; % d(yLinf)/dx
yL = yLinf;

% Integration: (Initialization overwritten for i > 2)
for i = 2:N-1
vL(i+1) = vL(i-1) + 2*dx*(exp(yL(i))-1);
yL(i+1) = yL(i-1) + 2*dx*vL(i);
end

% THE RIGHT HAND BOUNDARY CONDITION SOLUTION

% RHBC
A = -0.710 % Trial and error
yRinf = A*exp(-x); % RH asymtote
vRinf = -yRinf;

% Initialization
yR = yRinf;
vR = vRinf;

% Integration: (Initialization overwritten for i < N-1)
for i = N-1:-1:2
vR(i-1) = vR(i+1) - 2*dx*(exp(yR(i))-1);
yR(i-1) = yR(i+1) - 2*dx*vR(i);
end

% CONNECTING LHBC and RHBC SOLUTIONS

[dvymin imatchvy] = min(((vL-vR)*dx).^2 + (yL-yR).^2) %171

J = imatchvy
V = [vL(1:J-1); (vL(J)+vR(J))/2; vR(J+1:end)];
Y = [yL(1:J-1); (yL(J)+yR(J))/2; yR(J+1:end)];

k=k+1,figure(k)
subplot(211), hold on
plot(x,vL,'g--')
plot(x,vR,'r--')
plot(x,V)
plot(x(J),V(J),'*')
plot(x,zeros(N,1),'k')
legend('LHBG solution','RHBC solution','Merged Solution',3)
ylabel('V')
title('Solution to dV/dx = exp(Y) -1, dY/dx = V')
subplot(212), hold on
plot(x,yL,'g--')
plot(x,yR,'r--')
plot(x,Y)
plot(x(J),Y(J),'*')
plot(x,zeros(N,1),'k')
ylim([floor(min(Y))-1 ceil(max(Y)+1)])
legend('LHBG solution','RHBC solution','Merged Solution',4)
xlabel('x')
ylabel('Y')

Very Interesting Problem. Could probably find A
automatically.

Greg
First  |  Prev  | 
Pages: 1 2
Prev: How to get Fourier descriptor
Next: (Bode)2 figure in one