http://nl.mathworks.com/matlabcentral/newsreader/view_thread/317781
MATLAB Central Newsreader  ode, tan, singularity?
Feed for thread: ode, tan, singularity?
enus
©19942017 by MathWorks, Inc.
webmaster@mathworks.com
MATLAB Central Newsreader
http://blogs.law.harvard.edu/tech/rss
60
MathWorks
http://nl.mathworks.com/images/membrane_icon.gif

Thu, 08 Mar 2012 22:30:45 +0000
ode, tan, singularity?
http://nl.mathworks.com/matlabcentral/newsreader/view_thread/317781#869429
Roland
I want to plot a curve with a constant curvature. this code works very well:<br>
<br>
function curv()<br>
data.beta_2=89*pi/180;<br>
data.beta_1=20*pi/180;<br>
data.z_2 = 0;<br>
data.z_1 = 1;<br>
opts = odeset('events',@(u,z) cross_z_z_min(u,z,data));<br>
[u,z]=ode45(@(u,z) fun_dz_du(u,z,data),[0 inf],data.z_1,opts);<br>
plot(u,z)<br>
axis equal<br>
<br>
function dz_du = fun_dz_du(u,z,data)<br>
beta_2=data.beta_2;<br>
beta_1=data.beta_1;<br>
z_2 = data.z_2;<br>
z_1 = data.z_1;<br>
beta = beta_1  (beta_1beta_2)*(zz_1)/(z_2z_1);<br>
dz_du = tan(beta);<br>
<br>
function [gstop,isterminal,direction]=cross_z_z_min(u,z,data)<br>
gstop = zdata.z_2;<br>
isterminal = 1;<br>
direction = [];<br>
<br>
<br>
whereas when I change beta_2 to 91 degree the ode solver hangs up. I guess its related to the singularity of the tan function at 90 degree, but I dont know how to overcome that problem. Any ideas?

Fri, 09 Mar 2012 03:52:12 +0000
Re: ode, tan, singularity?
http://nl.mathworks.com/matlabcentral/newsreader/view_thread/317781#869451
Roger Stafford
"Roland " <burgmann@gmx.de> wrote in message <jjbbul$r7c$1@newscl01ah.mathworks.com>...<br>
> I want to plot a curve with a constant curvature. this code works very well:<br>
> <br>
> function curv()<br>
> data.beta_2=89*pi/180;<br>
> data.beta_1=20*pi/180;<br>
> data.z_2 = 0;<br>
> data.z_1 = 1;<br>
> opts = odeset('events',@(u,z) cross_z_z_min(u,z,data));<br>
> [u,z]=ode45(@(u,z) fun_dz_du(u,z,data),[0 inf],data.z_1,opts);<br>
> plot(u,z)<br>
> axis equal<br>
> <br>
> function dz_du = fun_dz_du(u,z,data)<br>
> beta_2=data.beta_2;<br>
> beta_1=data.beta_1;<br>
> z_2 = data.z_2;<br>
> z_1 = data.z_1;<br>
> beta = beta_1  (beta_1beta_2)*(zz_1)/(z_2z_1);<br>
> dz_du = tan(beta);<br>
> <br>
> function [gstop,isterminal,direction]=cross_z_z_min(u,z,data)<br>
> gstop = zdata.z_2;<br>
> isterminal = 1;<br>
> direction = [];<br>
> <br>
> whereas when I change beta_2 to 91 degree the ode solver hangs up. I guess its related to the singularity of the tan function at 90 degree, but I dont know how to overcome that problem. Any ideas?<br>
         <br>
This differential equation can be solved without using 'ode15'. The solution (if I understand your code correctly) is:<br>
<br>
z = 1/(b2b1)*(b2asin(sin(b1)*exp((b2b1)*u)))<br>
<br>
where b1 = data.beta_1 and b2 = data.beta_2, and ode45 is set to terminate when z descends to 0. If you have b2 > pi/2, then before the termination can occur the quantity inside the arcsine will eventually exceed 1 at which point it all blows up, but if b2 < pi/2 then ode45 terminates first before that can happen. Without your termination option, it would have hung up with either value of b2. The only way I can think of to "overcome" such a hang up is to change your criterion as to when ode45 is to stop. Stop it before b2(b2b1)*z gets past pi/2.<br>
<br>
What I don't understand is why you claim that this curve has a constant curvature. The only curves with constant curvature are circles and this curve is surely not a circle.<br>
<br>
Roger Stafford

Fri, 09 Mar 2012 09:03:36 +0000
Re: ode, tan, singularity?
http://nl.mathworks.com/matlabcentral/newsreader/view_thread/317781#869461
Roland
"Roger Stafford" wrote in message <jjbupc$m25$1@newscl01ah.mathworks.com>...<br>
> "Roland " <burgmann@gmx.de> wrote in message <jjbbul$r7c$1@newscl01ah.mathworks.com>...<br>
> > I want to plot a curve with a constant curvature. this code works very well:<br>
> > <br>
> > function curv()<br>
> > data.beta_2=89*pi/180;<br>
> > data.beta_1=20*pi/180;<br>
> > data.z_2 = 0;<br>
> > data.z_1 = 1;<br>
> > opts = odeset('events',@(u,z) cross_z_z_min(u,z,data));<br>
> > [u,z]=ode45(@(u,z) fun_dz_du(u,z,data),[0 inf],data.z_1,opts);<br>
> > plot(u,z)<br>
> > axis equal<br>
> > <br>
> > function dz_du = fun_dz_du(u,z,data)<br>
> > beta_2=data.beta_2;<br>
> > beta_1=data.beta_1;<br>
> > z_2 = data.z_2;<br>
> > z_1 = data.z_1;<br>
> > beta = beta_1  (beta_1beta_2)*(zz_1)/(z_2z_1);<br>
> > dz_du = tan(beta);<br>
> > <br>
> > function [gstop,isterminal,direction]=cross_z_z_min(u,z,data)<br>
> > gstop = zdata.z_2;<br>
> > isterminal = 1;<br>
> > direction = [];<br>
> > <br>
> > whereas when I change beta_2 to 91 degree the ode solver hangs up. I guess its related to the singularity of the tan function at 90 degree, but I dont know how to overcome that problem. Any ideas?<br>
>          <br>
> This differential equation can be solved without using 'ode15'. The solution (if I understand your code correctly) is:<br>
> <br>
> z = 1/(b2b1)*(b2asin(sin(b1)*exp((b2b1)*u)))<br>
> <br>
> where b1 = data.beta_1 and b2 = data.beta_2, and ode45 is set to terminate when z descends to 0. If you have b2 > pi/2, then before the termination can occur the quantity inside the arcsine will eventually exceed 1 at which point it all blows up, but if b2 < pi/2 then ode45 terminates first before that can happen. Without your termination option, it would have hung up with either value of b2. The only way I can think of to "overcome" such a hang up is to change your criterion as to when ode45 is to stop. Stop it before b2(b2b1)*z gets past pi/2.<br>
> <br>
> What I don't understand is why you claim that this curve has a constant curvature. The only curves with constant curvature are circles and this curve is surely not a circle.<br>
> <br>
> Roger Stafford<br>
<br>
<br>
<br>
thank you Roger for your reply. <br>
<br>
you are right, this does not create a curve with constant curvature. But what i want is a curve with a prescribed beta distribution along z. <br>
In this example, the distribution is linear, but it should be able to have any shape. thats why i prefere to use the ode solver. <br>
<br>
I am not shure if i understood correctly your suggestion. should i split up the function, one part for <br>
beta < pi /2 <br>
and one for <br>
beta > pi/2?

Fri, 09 Mar 2012 20:22:15 +0000
Re: ode, tan, singularity?
http://nl.mathworks.com/matlabcentral/newsreader/view_thread/317781#869520
Roger Stafford
"Roland " <burgmann@gmx.de> wrote in message <jjch18$ens$1@newscl01ah.mathworks.com>...<br>
> .........<br>
> I am not shure if i understood correctly your suggestion. should i split up the function, one part for <br>
> beta < pi /2 <br>
> and one for <br>
> beta > pi/2?<br>
       <br>
When (b2b1)/(z2z1) < 0, as is the case with your equations, the trajectory followed by a solution to those differential equations assumes an infinite dz/du slope at beta = pi/2 and u cannot continue to increase. The natural path for the trajectory to follow from there on would be along a vertical mirror image of the path with u now decreasing and z continuing to decrease as it finally asymptotically approaches the beta = pi level with u approaching minus infinity. However with ode45 committed to advancing u, this is not possible, so my recommendation is to set ode45 to stop right at beta = pi/2. As you have it set up, ode45 cannot continue from there without "hanging up" and producing nonsense.<br>
<br>
A similar situation holds for the entire family of trajectories with all the various possible initial conditions. Whenever beta is equal to pi/2 or differs from that by a multiple of pi  that is, whenever tan(beta) is infinite  the trajectories will necessarily reverse the direction of u as they cross the corresponding horizontal line. If (b2b1)/(z2z1) > 0 these trajectories are reversed in the u direction with asymptotes occurring as u approaches plus infinity.<br>
<br>
Roger Stafford

Sat, 10 Mar 2012 09:49:12 +0000
Re: ode, tan, singularity?
http://nl.mathworks.com/matlabcentral/newsreader/view_thread/317781#869553
Roland
"Roger Stafford" wrote in message <jjdopn$ptc$1@newscl01ah.mathworks.com>...<br>
> "Roland " <burgmann@gmx.de> wrote in message <jjch18$ens$1@newscl01ah.mathworks.com>...<br>
> > .........<br>
> > I am not shure if i understood correctly your suggestion. should i split up the function, one part for <br>
> > beta < pi /2 <br>
> > and one for <br>
> > beta > pi/2?<br>
>        <br>
> When (b2b1)/(z2z1) < 0, as is the case with your equations, the trajectory followed by a solution to those differential equations assumes an infinite dz/du slope at beta = pi/2 and u cannot continue to increase. The natural path for the trajectory to follow from there on would be along a vertical mirror image of the path with u now decreasing and z continuing to decrease as it finally asymptotically approaches the beta = pi level with u approaching minus infinity. However with ode45 committed to advancing u, this is not possible, so my recommendation is to set ode45 to stop right at beta = pi/2. As you have it set up, ode45 cannot continue from there without "hanging up" and producing nonsense.<br>
> <br>
> A similar situation holds for the entire family of trajectories with all the various possible initial conditions. Whenever beta is equal to pi/2 or differs from that by a multiple of pi  that is, whenever tan(beta) is infinite  the trajectories will necessarily reverse the direction of u as they cross the corresponding horizontal line. If (b2b1)/(z2z1) > 0 these trajectories are reversed in the u direction with asymptotes occurring as u approaches plus infinity.<br>
> <br>
> Roger Stafford<br>
<br>
<br>
<br>
Roger, thanks for the explanation of whats going wrong in my code. I solved the problem by solving for du/dz instead of dz/du, which makes everything much easier.<br>
<br>
Roland

Sat, 10 Mar 2012 14:35:12 +0000
Re: ode, tan, singularity?
http://nl.mathworks.com/matlabcentral/newsreader/view_thread/317781#869564
Roger Stafford
"Roland " <burgmann@gmx.de> wrote in message <jjf82o$g3m$1@newscl01ah.mathworks.com>...<br>
> Roger, thanks for the explanation of whats going wrong in my code. I solved the problem by solving for du/dz instead of dz/du, which makes everything much easier.<br>
         <br>
Yes, that would allow you to go past the beta = pi/2 point without stopping. Unfortunately, it will encounter severe accuracy problems in computing u when beta is near 0 or pi because cot(beta) approaches infinity there.<br>
<br>
Probably it would be preferable to allow arclength, s, to be the independent variable and send two simultaneous equations to ode45: either<br>
<br>
beta = .... % linear or otherwise in z<br>
du/ds = cos(beta)<br>
dz/ds = sin(beta)<br>
<br>
or else<br>
<br>
beta = ....<br>
du/ds = cos(beta)<br>
dz/ds = sin(beta)<br>
<br>
depending on which direction you wish s to be increasing. Either way you will still obtain the same curve determined by dz/du = tan(beta), and at the same time you have (du/ds)^2+(dz/ds)^2 = 1 showing that s is a genuine arclength along the curve. There should be no accuracy problems for beta either at pi/2 or near 0 and pi. An entire trajectory can be traced out within a [0,pi] range for beta without stopping (with the understanding that the curve is infinitely long in arclength.)<br>
<br>
Roger Stafford

Sat, 10 Mar 2012 16:49:21 +0000
Re: ode, tan, singularity?
http://nl.mathworks.com/matlabcentral/newsreader/view_thread/317781#869573
Roland
"Roger Stafford" wrote in message <jjfor0$3r6$1@newscl01ah.mathworks.com>...<br>
> "Roland " <burgmann@gmx.de> wrote in message <jjf82o$g3m$1@newscl01ah.mathworks.com>...<br>
> > Roger, thanks for the explanation of whats going wrong in my code. I solved the problem by solving for du/dz instead of dz/du, which makes everything much easier.<br>
>          <br>
> Yes, that would allow you to go past the beta = pi/2 point without stopping. Unfortunately, it will encounter severe accuracy problems in computing u when beta is near 0 or pi because cot(beta) approaches infinity there.<br>
> <br>
> Probably it would be preferable to allow arclength, s, to be the independent variable and send two simultaneous equations to ode45: either<br>
> <br>
> beta = .... % linear or otherwise in z<br>
> du/ds = cos(beta)<br>
> dz/ds = sin(beta)<br>
> <br>
> or else<br>
> <br>
> beta = ....<br>
> du/ds = cos(beta)<br>
> dz/ds = sin(beta)<br>
> <br>
> depending on which direction you wish s to be increasing. Either way you will still obtain the same curve determined by dz/du = tan(beta), and at the same time you have (du/ds)^2+(dz/ds)^2 = 1 showing that s is a genuine arclength along the curve. There should be no accuracy problems for beta either at pi/2 or near 0 and pi. An entire trajectory can be traced out within a [0,pi] range for beta without stopping (with the understanding that the curve is infinitely long in arclength.)<br>
> <br>
> Roger Stafford<br>
<br>
<br>
<br>
you are right, using arclength s as independent variable and calculating du/ds and dz/ds separately makes it more robust when beta approaches 0 and pi. I will implement it that way. <br>
Thanks a lot for you help!