30 views (last 30 days)

Hi everybody!

I've written a code that solves the following PB problem

This is a 2nd order ode and the dependant variable here is psi (tilde) and the independant variable is x (tilde). The boundary conditions are:

The code includes an optimizer based on the secant method. I'm trying to guess the solution that will ''land'' on the 2nd boundary condition via the shooting method:

function [x2,p2]=optimumup(y21,y22,bc2,l)

%optimizing the choice of y2 for the shooting method, using the Newton

%method (secant) by solving the equation p(length(p(:,1))-bc2=0

[x1,p1] = BP(y21,l);

[x2,p2] = BP(y22,l)

size(p2)

for i=1:10

dp=((p2(length(p2(:,1)),1)-bc2)-(p1(length(p1(:,1)),1)-bc2))/(y22-y21);

y23=y22-(p2(length(p2(:,1)),1)-bc2)/dp;%recurrence relation

p1=p2;

[x2,p2] = BP(y23,l);

if abs(p2(length(p2(:,1)))-bc2)<0.001;

disp(p2(length(p2(:,1)),1));

disp(p2(length(p2(:,1)),1)-bc2);

plot(x2,p2(:,1),x2,p2(:,2))

return

end

y21=p1(1,2);

y22=p2(1,2);

end

disp('No success')

end

function [x,p] = BP(y2,l)

%Shooting method for the BP eq

options = odeset('RelTol',1e-4,'AbsTol',1e-4);

[x,p] = ode45(@eqs,[0 l],[90/24 y2],options);

%n=size(p);p(n(1),1)

%plot(x,p(:,1),x,p(:,2))

end

function [dy] = eqs(x,y)

%PB equation

dy = zeros(2,1); % a column vector

dy(1) = y(2);

dy(2) = sinh(y(1));

end

y21 and y22 are the initial x guesses for the secant method. bc2 is the boundary condition at x goes to inifinity. l is the interval size.

I have two main problems:

- I don't know how to incorperate the inifinite boundary condition.
- For some reason my solution does not converge. I'm working dimensionless.

I would appreciate some help guys.

Thanks,

Roi

Alan Stevens
on 23 Oct 2020

I suspect the initial gradient required to get zero at infinity is an irrational number (or possibly a rational number with a large number of digits in its decimal expansion!). See the following for example:

xf = 25;

xspan = [0 xf];

dydx0 = -1.042119689394;

dx = 10^-11;

dydx0hi = dydx0+dx;

dydx0lo = dydx0-dx;

Y0 = [1, dydx0];

[x,Y] = ode45(@eqs,xspan,Y0);

psi = Y(:,1);

Y0hi = [1, dydx0hi];

[xhi,Yhi] = ode45(@eqshi,xspan,Y0hi);

psihi = Yhi(:,1);

Y0lo = [1, dydx0lo];

[xlo,Ylo] = ode45(@eqslo,xspan,Y0lo);

psilo = Ylo(:,1);

figure(1)

plot(x,psi,'k',xhi,psihi,'r',xlo,psilo,'b'),grid

axis([0 xf -2 2])

xlabel('x'),ylabel('\psi')

lbl1 = ['d\psi dx(x=0) = ',num2str(dydx0,12)];

lbl2 = ['d\psi dx(x=0) = ',num2str(dydx0hi,12)];

lbl3 = ['d\psi dx(x=0) = ',num2str(dydx0lo,12)];

legend(lbl1,lbl2,lbl3)

function dY = eqs(~,Y)

%PB equation

dY = [Y(2);

sinh(Y(1))];

end

function dY = eqshi(~,Y)

%PB equation

dY = [Y(2);

sinh(Y(1))];

end

function dY = eqslo(~,Y)

%PB equation

dY = [Y(2);

sinh(Y(1))];

end

This produces the following graph:

Although it might look as though the black line has converged on zero, if you extend the end time to, say 30, you will see it start to diverge.

Alan Stevens
on 24 Oct 2020

In effect I used the shooting method, except that I did it manually rather than writing a routine to do it! It’s actually not difficult, and took much less time than I would have taken to write a routine. Had the result looked as though the curve was tending to zero in a stable manner, I would then have written a routine to do it more accurately. However, as you can see from the fact that small perturbations in the initial slope lead to divergences in opposite directions, it seems that the final equilibrium is likely to be unstable.

A boundary condition at infinity suggests an asymptotic approach to zero (in this case). In physical systems, such as fluid flow boundary layers, for example, it is a convenient fiction standing in for a finite distance, so can often be replaced by a (relatively) large but finite value. If that is the case for you then you might be able to make use of the result. However, your equation doesn’t look as though it represents a physical system (I could be wrong!), so you will have to decide if it is useful or not.

Alan Stevens
on 24 Oct 2020

David Goodmanson
on 24 Oct 2020

Edited: David Goodmanson
on 24 Oct 2020

Hi Roi

For convenience I will use u in place of psi. An exact solution is

tanh(u/4) = A*exp(-x) where A = tanh(1/4) [1]

so

u = 4*atanh(A*exp(-x))

which meets both of the boundary conditons. For du/dx at the origin, it turns out that

du/dx = -2*sinh(1/2) = -1.042190610987495

which disagrees slightly from what Alan has, for reasons not known to me.

===========derivation=============

start with

u'' = sinh(u)

multiply by integrating factor u' and integrate

u'^2/2 = cosh(u) (+ constant)

u'^2 = 2*cosh(u) + constant

as x --> inf you want u --> 0 which also means that u' --> 0, so the constant = -2

u'^2 = 2*(cosh(u)-1) = 4*sinh(u/2)^2

u' = -2*sinh(u/2) % pick the negatatve root

du/(2*sinh(u/2)) = -dx

% integrate

log(tanh(u/4)) = -x + constant

tanh(u/4) = A*exp(-x)

you want u = 1 at x = 0 so

A = tanh(1/4)

For the derivative at the origin, take the derivative of both sides of [1] at the top

( (1/4) / cosh^2(u/4) ) du/dx = -A*exp(-x)

% evaluate this at x = 0 and u(0) = 1

( (1/4) / cosh^2(1/4) ) du/dx = -A

du/dx = -tanh(1/4)*4*cosh^2(1/4) = -4*sinh(1/4)*cosh(1/4) = -2*sinh(1/2)

David Goodmanson
on 24 Oct 2020

Hi Alan,

Thanks, it did work in this special case but I think your comments about a finite boundary and/or an asymptotic approach are going to be more useful in general. In this case out at infinity the answer has to be u = A*exp(-x), so both u and u' are known there. Then starting at large x, integrating to the origin and using the shooting method for various A should also work, I think. I have not tried it.

Roi Bar-On
on 25 Oct 2020 at 15:16

Alan Stevens
on 30 Oct 2020 at 22:32

Opportunities for recent engineering grads.

Apply TodayFind the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!
## 0 Comments

Sign in to comment.