How to restrict the domain of dependent variable in ode45

Hello every one

I have this equation of ray evolution that I want to solve in matlab using ode45

where dn2 / dx equals to n * (dn / dx) and n is a function of x as n(x). n is the refractive index in Optics and this equation shows how the ray evolves in a medium of height x and range z.

β is a constant and existing, but here is my problem: I do not have the explicit form of n(x), instead, I have a data set of n for x = (0,30). How can I solve this equation using ode45, or any other method compatible for this problem, in a way that the dependent variable x is considered between 0 to 30, because my n varibales are bound only to this interval, so that I would only substitute the n variables and rest assured that x is changing according to n. n is a 1 * 300 matrix by the way

[J. Alex Lee]

are you sure about the equation itself? what aspect of the "ray" is it supposed to describe? I'm not familiar with optics, but it seems odd that you are describing x and z as space coordinates, and yet your ode is of the form x''[z] = ...

[Proman]

Yes sir, the equation is fully correct. It describes the ray path in a medium

[Proman]

the order of the equation does not matter since it is reducible to a system of first order ode as:

function dx = odefun(z,x)

dx(1) = x(2);

dx(2) = 1/2β^2 * n * (dn / dx)

end

[James Tursa]

If n is only a data set, then it seems your option will be to curve fit and then interpolate to come up with the n values. Differentiate the curve fit and then interpolate for the dn/dx values.

[J. Alex Lee]

agree with James Tursa...what does your n(x) data look like? The core of your problem appears to be that you don't know what n(x) is outside the range 0<x<30, but you really do need that unless you are really lucky and the combination of BCs keeps you in the range. Depending on your your n(x) data looks, maybe you can make reasonable extrapolations, and just check a posteriori what is the range of your solution.

[Proman]

Well, I reduced my problem to the following equation; All I need to do is numerically solve this equation; how can I numerically solve it?

[Walter Roberson]

Is your still an unknown function that needs to be deduced from a limited set of data?

I am not clear from what you were saying before that you know the boundary conditions n(0) and n(30), or if you were saying that you know n(0), n(1), ... n(30), or if you were saying that you know n at some vector of values that are all in the range 0 to 30 ?

And is it correct that x is a function of z, so n(x) is really n(x(z)) ? Is there reason to believe that x(z) passes through the area over which you know some n values?

[J. Alex Lee]

If your reduced equation is correct, it seems to me at first glance to imply that you actually know more about the interrelations of n, x, z than you stated...otherwise I don't think that is generally equivalent to the first equation...

Walter's first part of last question is what I was asking before; it just seems an odd equation, but I admit ignorance in the application area.

[Proman]

I attached the file here. n is a 1*300 vector obtained from a function in my code called BYCModelIRfun; In this function, z is from 0.1 to 30 so n(z) is for these heights. Then I solve this equation for different n but I do not know why I do not get the correct answer.

[J. Alex Lee]

this doesn't help..why not just give the resulting data. is "n(z)" here a typo for "n(x)", or the other way around? that would be a pretty important detail

[Proman]

yes n(z) is n(x). n is a function of the dependent variable x which is n(x). that was my mistake. I attached the .mat data (n.mat)

the initial condition is x(0) = 10. n(x) data are for x = (0.1,30). z span is z = (0,10000).

suppose β (denoted by "B") equals to 1.000256103908571. this is my proposed algorithm which I believe is wrong:

x(1,1) = 10;
B = 1.000256103908571;
z = linspace(0,10000,length(n));
for s = 1 : length(n) - 1
dz = z(s + 1) - z(s);
dx(1,s) = (sqrt((n(s) ^ 2 - B ^ 2) - 1)) * dz;
x(1,s+1) = dx(1,s) + x(1,s);
end

Thanks again for supporting me and helping me in this matter. Much appreciated Sir!

[J. Alex Lee]

you've confused me with the example code that you have posted...it doesn't appear to rely on ode45 at all...

but as interpolation/extrapolating, you might say with confidence that the function n(x) is linear beyond x = 30, and safely extrapolate-able. Not sure about the left side though. Can you compute for x<0.1 with your above code?

Alternatively, since you have the code that generates n values, can you package that into a function of the form

function n = n_func(x)
% ...
end

and just call it in your odefun?

[Proman]

n_func is the function generates n and Main reads n and trace the ray and plots it

I want to see if my algorithm is right or not

Note: In this code, the dependent variable is z and the independent one is x (unlike what is defined in the above equation).

Again thank you for your kind help

[J. Alex Lee]

The algorithm you posted above, which looks somewhat like a manual discretization? If the question is whether that algorithm gets you the solution to the ODE you posted, then I would say not at all. It looks like you are attempting forward Euler (why? why not ode45 as you suggested in original post), but it is nonsensical to me that you are defining your z domain based on the "length of the vector n".

[Proman]

Yes, the most simple algorithm for a 1st order ode is the one I did as manual discretiation. Well, I am not surprised that the algorithm proved wrong in my problem but how can I fix it to my interest Can you help me at it? I did not know how to use ode45 in this problem so I came up with that algorithm.

I defined z domain based on the length of vector n because I thought I had to plot z in terms of x so I had to discrrtize it accordingly but as you said, I am extremely interested that its steps does not necessarily equal to the size of n but I do not how.

the whole problem is eating me up. How can I get rid of this riddle? can you help me edit my Main.m code? I just want something to get me the rigth solution, no matter what method I am supposed to use; ode45 or any other method at hand.

[Walter Roberson]

I still do not understand what the question actually is. I do not see that my questions of https://www.mathworks.com/matlabcentral/answers/573340-how-to-restrict-the-domain-of-dependent-variable-in-ode45#comment_957058 were answered.

[Proman]

I answered your question in the commnets but maybe it was vague. I am sorry for my mistake. I have already attached my function (called n_func.m generating n(x)) and Main.m (calling n_func.m) which solves the equation. (I attached them again on this comment).

Let me explain from the scratch (and according to my code):

This is my equation:

here, z is the dependent variable (height) and x is the independent one (range). the purpose of the code is to simulate the ray path in the vicinty of the sea surface (for each point at range x and height z and for different angle of incident light which in the code are shown as i1); In optics, it is not traditionally convincing to say, for instance, z is a function of x and I completely understand what you mean because it is really not, but the equation implies that for tracing the ray, we have to consider such situation as if z is a function of x.

we do not have the explicit form of n(z) as a function (no polynimialor transcedental form); instead, we are going to obtain it numerically from the function n_func.m for z = (0.1,30). It means that, for every z, there will be a n(z). in this function, we discretized z to 300 points, so n(z) would be 300 points as well and n is n(z(x)). we can not set z = 0 becuase we have some logarithmic restriction, so we have to start from z = 0.1 or eps. Having n(z) at hand, we want to solve the equation above for x = (0,10000), and we want to start ray tracing at z(0) = 10 (which is our initial condition in the equation).

Considering these, How can I solve this equation to get the right answer? I already proposed a very naive and simple algorithm which it seems to be wrong.

You also asked "Is there reason to believe that z(x) passes through the area over which you know some n values? ". that is my quaestion too to be honest. In fact I have no opinion if your question is necessary to solve this problem or not but I believe that your doubt might be right and we maybe need to extrapolate n, yet I insist that I am not sure.

Thank you for your attention so far and I am very happy that you followed up until now. If anyother information needed, I am at your service.

[J. Alex Lee]

You have now flipped x's and z's so many times it's hard to know what's what. It seems Walter's question is the most pertinent.

If you have no reason to believe your solution will be bound to the range (0.1,30), then I would say simply test it a posteriori. Follow James's advice to interpolate, and see if your algorithm fails because you try to evaluate n(z) at some value outside the range.

At least for z>30, I've observed that you can probably safely extrapolate linearly.

If you know the asymptotic form of as , you could piece-wise define n(z) as an interpolant from , analytic form for , and linear from .

If β value is less than the trough-like feature in you will likely have no issues, as this will guarantee your derivative to be positive so that .

[Proman]

Then I am all set! Many thanks for your help ^-^