# Solving for Coefficient in Nonlinear 2nd-Order Differential Equation

1 view (last 30 days)
Alexander Herrod on 4 Jan 2018
Commented: Alexander Herrod on 22 Jul 2018
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!

Torsten on 5 Jan 2018
Edited: Torsten on 5 Jan 2018
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 on 21 Jul 2018
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!)

John BG on 5 Jan 2018
Edited: John BG on 15 Jan 2018
Hi Alexander
1.
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;
thanks in advance for time and attention
John BG
Alexander Herrod on 22 Jul 2018
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!