Soling a boundary value problem

I am trying to solve the following elliptic equation

on the interval [0,1] with Neumann Boundary condition, I want to get a positive solution φ with the maximum of φ is 1.

The following is the code I am using

h=0.001;

L=1;

M=L/h;

x=linspace(0,L,M+1);

gamma = @(x)(1+0.1*sin(2*pi*x)).*(1+0.1*sin(2*pi*x));

beta=@(x)(1+0.1*sin(2*pi*x));

%v0 = @(x)(1+0.5*cos(pi*x));

%v0x=@(x)-0.5*pi*sin(pi*x);

bvpfcn = @(x,v)[v(2);(gamma(x)-(beta(x)).*v(1)];

bcfcn = @(va,vb)[va(2);vb(2)];

guess = @(x)[1;0];

%guess = @(x)[v0(x);v0x(x)];

solinit = bvpinit(x,guess);

sol = bvp4c(bvpfcn, bcfcn, solinit);

phi1=sol.y(1,:);

M=max(phi1);

phi=phi1/M;

plot(sol.x,phi)

I have difficulity in selecting the initial function. I tried [1;0] and [v0(x);v0x(x)] which is defined above, I got different solutions. I cannot figure out the problem. Any help would be appreciated!

Answer 1

[Torsten]

Prove the following:

If u is a solution of your differential equation, so is c*u for every c in IR.

So it does not surprise that bvp4c gives different solutions dependent on your initial guess: the solution is not unique.

If you normalize the solution as you did after solving, it seems it becomes unique. You can see this if you solve your equation for your two initial guess functions and plot the solutions in one graph.

[Sharon]

I poltted them in one graph. They are not the same solution.

[Torsten]

It seems phi = 0 is the only solution you get for your equation. The differences in the solutions for different initial functions are in the order of 1e-12, as you can see in the plot below.

h=0.001;
L=1;
M=L/h;
x=linspace(0,L,M+1);
gamma = @(x)(1+0.1*sin(2*pi*x)).*(1+0.1*sin(2*pi*x));
beta=@(x)(1+0.1*sin(2*pi*x));
v0 = @(x)(1+0.5*cos(pi*x));
v0x=@(x)-0.5*pi*sin(pi*x);
bvpfcn = @(x,v)[v(2);(gamma(x)-beta(x))*v(1)];
bcfcn = @(va,vb)[va(2);vb(2)];
guess = @(x)[1;0];
solinit = bvpinit(x,guess);
options = odeset('RelTol',1e-8,'AbsTol',1e-8);
sol = bvp4c(bvpfcn, bcfcn, solinit,options);
hold on
plot(sol.x,sol.y(1,:))
guess = @(x)[v0(x);v0x(x)];
solinit = bvpinit(x,guess);
sol = bvp4c(bvpfcn, bcfcn, solinit,options);
plot(sol.x,sol.y(1,:))
hold off

[Sharon]

Thanks for your reply! phi=0 is a solution. I think there are positive solutions. I want to get a positive solution with the maximum 1. As you said before, if I normalize the solution, the solution is unique.

[Torsten]

How do you know there are solutions different from phi = 0 ? If one such solution phi0 would exist, there would be infinitly many of them because - as noticed above - also c*phi0 would be a solution.

[Sharon]

Yes, there are infinitely many solutions without normalization.

[Sharon]

Are there other methods to solve such ode?

[Torsten]

It seems that bvp4c as a numerical solver always converges to the trivial solution phi = 0.

I tried "dsolve", but it is unable to find an analytical solution (s.b.).

I'm still interested in your source that the equation has a solution different from 0.

syms x y(x)
eqn = diff(y,x,2)-((1+0.1*sin(2*pi*x))^2-(1+0.1*sin(2*pi*x)))*y==0

eqn(x) = 

dsolve(eqn)

Warning: Unable to find symbolic solution.

 
ans =
 
[ empty sym ]
 

[Sharon]

Thank you very much for your time and efforts!