ode45 does not solve my second order problem

Hello to everyone

I have a second order ode that I want to solve for this equation (I need to say that I already attached my code to this post in the last section):

where β^2 is a constant and

where n0, n2 and alpha are constant values and all of these constants can be found in my code attached here. when I try to solve this ode equation with ode45, I get a warning as "Warning: Failure at t=3.069827e+00. Unable to meet integration tolerances without reducing..." and my answer in the plot does not show the one I need to get. this is what I want to get:

and this is my result:

which you could make sure by running my code, If you wanted to. how can I solve this problem. I would appreciate any help and thanks in advance for your time. I need to say that for other forms of n(z), matlab exactly solved the equation as I wanted but this one is tricky and I can not figure out in what way. here is my code by the way if you can not download the file

function example2
clc
clear 
for i1 = [pi/10 0 -pi/60 -pi/30 -pi/15 -pi/10]
xspan = [0 20] ;
z0=[0.5,tan(i1)] ;
hold on

[x,z] = ode45(@odefun,xspan,z0);
plot(x,z(:,1))
xlabel('Range(m)')
ylabel('height(m)')
grid on
axis([0 20 0 2])

end

function dz=odefun(~,z)

a = 2.303;
    
n2 = 0.45836;

n0 = 1.000233;

np = sqrt(n0 ^ 2 + (n2 ^ 2 * exp(-a * 0.5)));
 
B = np * cos(i1); 

n = sqrt(n0 ^ 2 + (n2 ^ 2 * exp(-a * z(1))));

dn_dz = -(a*n2^2*exp(-a*z(1)))/(2*(n0^2 + exp(-a*z(1))*n2^2)^(1/2));

dz=zeros(2,1); 

dz(1)=z(2);

dz(2)= (n) * (dn_dz) / (B^2);


end

end

Answer 1

[Alan Stevens]

Your warnings occur when z(1) goes negative. They can be eliminated by including

if z(1)<0
   z(1) = 0; 
end

at the beginning of your odefun function.

Are your sure your initial angles are the same as the plot from the text? Try using these angles instead

i1 =[pi/13 pi/15 pi/20 0 -pi/60 -pi/30 -pi/15 -pi/10]

You should get

[Proman]

Many thanks my friend. Absolutely correct.

I think you are right and the angles need to be modified.

Answer 2

[Walter Roberson]

This has exactly the same problem as before. Your left hand side involves z as a function of x, but your right hand side involves z as an independent variable, not as a function of x.

You even give a definition for n(z) that does not involve x in any obvious way.

If we say "Okay, so z is really z(x)" then the is a derivative with respect to a function, which is not valid in normal calculus (requires Calculus of Variations.)

In order to be able to differentiate n with respect to z even though z is a function of x, then you would need to know that n would be considered independent of z for calculus purposes, which is something you cannot know because z(x) is not known.

However... in theory you can try it, as long as you back-substitute the z(x) you find, to check to see that the derivative you get is consistent with the hypothesized derivative.

syms beta n n_0 n_2 alpha z(x) Z
n = sqrt(n_0^2 + n_2^2 * exp(-alpha*z(x)));
dndx = subs(diff(subs(n,z,Z),Z),Z,z) * diff(z(x),x);
disp(n)
disp(dndx)
d2zdx = diff(z(x),x,x);
disp(d2zdx)
eqn = d2zdx == n / beta^2 * dndx;
disp(eqn);
ZX = dsolve(eqn)

ZX =
 log((- n_2^2 + exp(C1*alpha*(C2 + x)))/(2*C1*beta^2))/alpha
                                                          C1
 

Two solutions for z(x), one of which is constant and the other is more interesting.

Now, can we use the boundary conditions from your code in order to proceed to find the values of C1 and C2 based upon the boundary conditions?

Well.. NO. Or at least not easily. You have the boundary condition for n(0) = 1/2 . In order to phrase that in terms of x so that you can found boundary conditions for z(x), you would have to solve z(x) == 0 to find x

Z0 = solve(subs(n,z,ZX(1))==0,x);   %this will warn about parameterizing by k
vars = symvar(ZX);   %C1, C2, alpha, beta, n_2, x
C2 = solve(subs(n,z,Z0) == 1/2,vars(2));
dZ0 = simplify(subs(subs(subs(dndx, z, ZX(1)),vars(2),C2),x,0));
C1 = solve(dZ0==tan(Pi/10),vars(1))

and there we run out of steam: the equation is too messy to find the C1 constant.

"where β^2 is a constant"

np = sqrt(n0 ^ 2 + (n2 ^ 2 * exp(-a * 0.5)));
B = np * cos(i1); 

That doesn't look like a constant to me.

[Proman]

Thank you for your time and follow-up on my question. It is really kind of you. To be honest, you are quite a master on these things! However, in my problem, n is not a function of x. In my previous question, n was a data but here it is a function.

I exactly know what you try to point out yet things are not that complicated in optics, yet mathematically, you are srongly correct.

To be honest, It is a task I have to see it done and personally, I still have not make heads and tails of this equation and have more questions on it even more than you.

If you find this material interesting, I refer you to the exempt from a book I attached to this comment in "Chapter 3, Geometrical Optics" (specifically page 39, section 3.4: The ray equation and its solutions) and since you are a lot stronger than me in such mathematical material, you may better understand it.

If you had any questions regarding the optical background, It would be an honor to be at your service sir.