ode45 does not solve my second order problem

2 次查看（过去 30 天）
Proman
Proman 2020-8-9
评论： Proman 2020-8-9
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

请先登录，再回答此问题。

采纳的回答

Alan Stevens
Alan Stevens 2020-8-9
编辑：Alan Stevens 2020-8-9
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
  1 个评论
Proman
Proman 2020-8-9
Many thanks my friend. Absolutely correct.
I think you are right and the angles need to be modified.

请先登录，再进行评论。

更多回答（1 个）

Walter Roberson
Walter Roberson 2020-8-9
编辑：Walter Roberson 2020-8-9
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.
  1 个评论
Proman
Proman 2020-8-9
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.

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Ordinary Differential Equations 的更多信息

标签

产品


版本

R2019b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by