Solving for Coefficient in Nonlinear 2nd-Order Differential Equation

Hi,

I'm trying to numerically solve for the coefficient "a" in a nonlinear second-order differential equation of the form:

y''(r) = y(r)(a - b*abs(y(r))/r + (c*r^3 - d)*y(r)^2)

given boundary conditions:

y(0) = y(inf) = 0

integral from 0 to inf of y(r)^2 dr = 1

However I don't seem to be able to find the tools to do this neatly. I have another platform that can do this, but I'd very much like to cross-check the results in Matlab!

Any help is appreciated!

Answer 1

[Torsten]

Use "fzero" or "fsolve" to solve the equation

(integral_{0}^{inf} y(a,r)^2 dr) - 1 = 0

for a.

To this end, in the function routine from "fzero" or "fsolve", you will be given a suggestion for a. With this value for a, you can use "bvp4c" to solve the boundary value problem

y''(r) = y(r)(a - b*abs(y(r))/r + (c*r^3 - d)*y(r)^2)

given boundary conditions:

y(0) = y(inf) = 0

Once you have its solution in discrete points, use "trapz" to approximate integral_{0}^{inf} y(r)^2 dr and return integral_{0}^{inf} y(r)^2 dr - 1 to "fzero" or "fsolve".

I think one problem with a numerical solution will be that your boundary value problem always has y = 0 for all r as solution.

Best wishes

Torsten.

[Alexander Herrod]

Thanks for the reply, and sorry it's taken so long!

I don't assume the dependency of y(r) on any of the constants. The integration requirement is only to be applied at the end.

If I try your suggestion, Matlab throws me an error that it can't find and solutions to that boundary condition (because of course there are infinite possible functions!)

Answer 2

[John BG]

Hi Alexander

this is John BG <mailto:[e-mail] [e-mail]>

1.

One of the tools available is MuPAD, start MuPAD with

mupad

once open define the equation, but avoid the abs(), it only complicates the result and since r seems to be a cylindrical coordinate, if r is the radius, then there is no point that r takes negative values

.

the result after ode:solve(o) is what comes up if including the abs()

the resulting functions do not seem that useful.

2.

It turns out that the Jacobi elliptic functions are solution to a similar differential equation:

clear all;clc;close all

rspan=[-20 20];
dr=0.0001;
r=[rspan(1):dr:rspan(2)];

M = 0.7;
[S,C,D] = ellipj(r,M);
figure(1);plot(r,S,r,C,r,D);
legend('SN','CN','DN','Location','best');grid on;
title('Jacobi Elliptic Functions sn,cn,dn');

a=1;b=1;c=.1;d=.1;

vS=d2ydr2(r,S,a,b,c,d);
vC=d2ydr2(r,C,a,b,c,d);
vD=d2ydr2(r,D,a,b,c,d);

d2Sdt2=diff(diff(S));
d2Cdt2=diff(diff(C));
d2Ddt2=diff(diff(D));

errS=d2Sdt2-vS([3:end]);
errC=d2Cdt2-vC([3:end]);
errD=d2Ddt2-vD([3:end]);

figure(2);plot(r([3:end]),abs(errS),r([3:end]),abs(errC),r([3:end]),abs(errD));grid on;

3.

now M=.99 CN and DN (not the blue one) show far smaller error

rspan=[1 5];
dr=0.0001;
r=[rspan(1):dr:rspan(2)];
M = 0.99;
[S,C,D] = ellipj(r,M);

vS=d2ydr2(r,S,a,b,c,d);
vC=d2ydr2(r,C,a,b,c,d);
vD=d2ydr2(r,D,a,b,c,d);

d2Sdt2=diff(diff(S));
d2Cdt2=diff(diff(C));
d2Ddt2=diff(diff(D));

errS=d2Sdt2-vS([3:end]);
errC=d2Cdt2-vC([3:end]);
errD=d2Ddt2-vD([3:end]);
figure(3);plot(r([3:end]),abs(errS),r([3:end]),abs(errC),r([3:end]),abs(errD));grid on;

if you find this answer useful would you please be so kind to consider marking my answer as Accepted Answer?

To any other reader, if you find this answer useful please consider clicking on the thumbs-up vote link

thanks in advance for time and attention

John BG

<mailto:[e-mail] [e-mail]>

[Alexander Herrod]

Thanks John,

Sorry it's taken me so long to get to this. Is there any way to input boundary conditions of the kind I've given?

Thanks again!