Solve second order non linear differential equation

I need to solve F'' + F^2 -1/2pi = 0

Boundary conditons

F(0) = 0;

F(inf) = 1

I am new to using the ode solver in matlab and am not sure how to make it solve a non-linear SECOND order equation. Any suggestion would be appreciated.

Answer 1

[Sulaymon Eshkabilov]

You can try a builtin fcn: bvp4c(). Here is a nice documentation on bvp4c() with examples: https://www.mathworks.com/help/matlab/ref/bvp4c.html

[AJ]

Thanks!

[Sulaymon Eshkabilov]

Most Welcome!

[John D'Errico]

I don't think the bvpsolver can handle an infinite domain, since it uses collocation.

Answer 2

[John D'Errico]

You will not be able to use bvp4c directly, since it will use collocation. That will fail on an infinite domain.

One option might be to use a finite domain, but one that is relatively large. So [0,Xmax] where Xmax is sufficiently far out that the solve has a chance to finish in a reasonable time, yet not that far out that numerical problems exist.

I did try dsolve, and it fails to find a solution, but that is not unexpected. My next suggestion is to use a transformation of the domain into a bounded domain. For example, you might try a transformation like

y = atan(x)

This maps the domain [0,inf] into [0,pi/2], which may be surprisingly a good idea, considering the equation at hand. You can see some ideas here on a similar problem:

https://nicoguaro.github.io/posts/infinite_fdm/

But certainly other links exist where people explain how to transform such a problem into a finite domain.

[AJ]

Thanks John! I think I can use Xmax for a large value.

[John D'Errico]

As a quick shot, it might look like this.

xmax = 100;
xmesh = linspace(0,xmax,1001);
guess = @(x) [x/xmax;repmat(1/xmax,size(x))]; % linear, but satisfies the BC
solinit = bvpinit(xmesh, guess)

solinit = struct with fields:
    solver: 'bvpinit'
         x: [0 0.1000 0.2000 0.3000 0.4000 0.5000 0.6000 0.7000 0.8000 0.9000 1 1.1000 1.2000 1.3000 1.4000 1.5000 1.6000 1.7000 1.8000 1.9000 2 2.1000 2.2000 2.3000 2.4000 2.5000 2.6000 2.7000 2.8000 2.9000 3 3.1000 3.2000 3.3000 3.4000 3.5000 … ]
         y: [2×1001 double]
     yinit: @(x)[x/xmax;repmat(1/xmax,size(x))]

bcfun = @(ya,yb) [ya(1);yb(1) - 1]; % zero at x=zero, 1 at x=xmax

% the second order BVP, sonverted to a system of two first order equations.
% Since the ODE is:
%  y'' + y^2 - pi/2 = 0
%
% converted to a 1st order system, we will have
%  y1 = y
%  y2 = y'
%
% Then the problem reduces to the nonlinear system:
%  y1' = y2
%  y2' = -y1^2 + pi/2
bvpfun = @(X,Y) [Y(2);-Y(1)^2 + pi/2];

sol = bvp4c(bvpfun,bcfun,solinit)

Warning: Unable to meet the tolerance without using more than 5000 mesh points.
The last mesh of 2890 points and the solution are available in the output argument.
The maximum residual is 531.443, while requested accuracy is 0.001.

sol = struct with fields:
    solver: 'bvp4c'
         x: [0 0.0333 0.0667 0.1000 0.1500 0.2000 0.2333 0.2667 0.3000 0.3333 0.3667 0.4000 0.4333 0.4667 0.5000 0.5333 0.5667 0.6000 0.6333 0.6667 0.7000 0.7333 0.7667 0.8000 0.8333 0.8667 0.9000 0.9333 0.9667 1 1.0333 1.0667 1.1000 1.1333 1.1667 … ]
         y: [2×2890 double]
        yp: [2×2890 double]
     stats: [1×1 struct]

plot(sol.x,sol.y)

Which does not look very satisfying to me. I cannot really talk about any limiting behavior as x --> inf, since it seems to be oscillating between two limits for large x. And that means what it returned as a possible solution is pretty much pure crapola. I do note that bvp4c did not think it had converged. That might suggest I got something wrong in the equations (easy enough, I was not terribly careful) or that I provided a poor initial guess.

So perhaps a better guess might make sense. Suppose I guess the relationship looks like sin(x) over that range, thus starts at zero, then approaches 1 at the right end point?

guess = @(x) [sin(x/xmax*pi/2);cos(x/xmax*pi/2)*pi/2/xmax];
solinit = bvpinit(xmesh, guess);
sol = bvp4c(bvpfun,bcfun,solinit);

Warning: Unable to meet the tolerance without using more than 5000 mesh points.
The last mesh of 2179 points and the solution are available in the output argument.
The maximum residual is 164.824, while requested accuracy is 0.001.

plot(sol.x,sol.y)

Far less nasty looking, but still not happy.

Ok, consider if my initial guess looks like erf? Erf has the necessary properties, that it satisfies the bounday conditions, and it is trivial to differentiate.

xmax = 20;
xmesh = linspace(0,xmax,1001);
guess = @(x) [erf(x);1/sqrt(pi)*exp(-x.^2)]; % linear, but satisfies the BC
solinit = bvpinit(xmesh, guess);
plot(solinit.x,solinit.y)

sol = bvp4c(bvpfun,bcfun,solinit);

At least this time, bvp4c does not have a hissy fit about convergence. Probably due to the narrower interval.

plot(sol.x,sol.y)

Again, while I chose a smaller interval, the solution bvp4c finds is a purely oscillatory one. And that makes the boundary condition as x--> inf would seem to be meaningless if the function is a sinusoidally oscillating one. Actually though, it is not so, since you can have a function that is everywhere oscillatory, yet has a finite limit at infinity. For example:

y = 1 + exp(-x)*sin(x)

has a limit of 1 as x--> inf, yet is purely oscillatory. Is that the behavior I would expect for any solution of this ODE? Well, actually, yes. That is, suppose F(x) is a function that approaches a limit of 1 at inf, where the derivatives, so F'(x) and F''(x) all approach zero? That cannot be a solution to your ODE, thus:

F''(x) + F(x)^2 - pi/2 = 0

since then as x--> inf, if F''(x)--0, then we would have

1^2 - pi/2 = 0

which is clearly false. And that means ANY solution to this problem is indeed an oscillatory one. That relieves me, since bvp4c was trying to tell me that all along.

[AJ]

Thank you!John! Your efforts are appreciated

Answer 3

[MOSLI KARIM]

Pi=1;
BC=@(ya,yb)[ya(1)
    yb(1)-1];
dxdy=@(x,y)[y(2);(1/2)*Pi-y(1)^2];
solinit=bvpinit(linspace(0,1,5),[0 0]);
sol=bvp4c(dxdy,BC,solinit);
figure(2)
hold on 
plot(sol.x,sol.y)