Discover MakerZone

MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

Learn more

Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

Apply Today

Thread Subject:
ode, tan, singularity?

Subject: ode, tan, singularity?

From: Roland

Date: 8 Mar, 2012 22:30:45

Message: 1 of 7

I want to plot a curve with a constant curvature. this code works very well:

function curv()
data.beta_2=89*pi/180;
data.beta_1=20*pi/180;
data.z_2 = 0;
data.z_1 = 1;
opts = odeset('events',@(u,z) cross_z_z_min(u,z,data));
[u,z]=ode45(@(u,z) fun_dz_du(u,z,data),[0 inf],data.z_1,opts);
plot(u,z)
axis equal

function dz_du = fun_dz_du(u,z,data)
beta_2=data.beta_2;
beta_1=data.beta_1;
z_2 = data.z_2;
z_1 = data.z_1;
beta = beta_1 - (beta_1-beta_2)*(z-z_1)/(z_2-z_1);
dz_du = -tan(beta);

function [gstop,isterminal,direction]=cross_z_z_min(u,z,data)
gstop = z-data.z_2;
isterminal = 1;
direction = [];


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?

Subject: ode, tan, singularity?

From: Roger Stafford

Date: 9 Mar, 2012 03:52:12

Message: 2 of 7

"Roland " <burgmann@gmx.de> wrote in message <jjbbul$r7c$1@newscl01ah.mathworks.com>...
> I want to plot a curve with a constant curvature. this code works very well:
>
> function curv()
> data.beta_2=89*pi/180;
> data.beta_1=20*pi/180;
> data.z_2 = 0;
> data.z_1 = 1;
> opts = odeset('events',@(u,z) cross_z_z_min(u,z,data));
> [u,z]=ode45(@(u,z) fun_dz_du(u,z,data),[0 inf],data.z_1,opts);
> plot(u,z)
> axis equal
>
> function dz_du = fun_dz_du(u,z,data)
> beta_2=data.beta_2;
> beta_1=data.beta_1;
> z_2 = data.z_2;
> z_1 = data.z_1;
> beta = beta_1 - (beta_1-beta_2)*(z-z_1)/(z_2-z_1);
> dz_du = -tan(beta);
>
> function [gstop,isterminal,direction]=cross_z_z_min(u,z,data)
> gstop = z-data.z_2;
> isterminal = 1;
> direction = [];
>
> 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?
- - - - - - - - - -
  This differential equation can be solved without using 'ode15'. The solution (if I understand your code correctly) is:

 z = 1/(b2-b1)*(b2-asin(sin(b1)*exp((b2-b1)*u)))

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-(b2-b1)*z gets past pi/2.

  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.

Roger Stafford

Subject: ode, tan, singularity?

From: Roland

Date: 9 Mar, 2012 09:03:36

Message: 3 of 7

"Roger Stafford" wrote in message <jjbupc$m25$1@newscl01ah.mathworks.com>...
> "Roland " <burgmann@gmx.de> wrote in message <jjbbul$r7c$1@newscl01ah.mathworks.com>...
> > I want to plot a curve with a constant curvature. this code works very well:
> >
> > function curv()
> > data.beta_2=89*pi/180;
> > data.beta_1=20*pi/180;
> > data.z_2 = 0;
> > data.z_1 = 1;
> > opts = odeset('events',@(u,z) cross_z_z_min(u,z,data));
> > [u,z]=ode45(@(u,z) fun_dz_du(u,z,data),[0 inf],data.z_1,opts);
> > plot(u,z)
> > axis equal
> >
> > function dz_du = fun_dz_du(u,z,data)
> > beta_2=data.beta_2;
> > beta_1=data.beta_1;
> > z_2 = data.z_2;
> > z_1 = data.z_1;
> > beta = beta_1 - (beta_1-beta_2)*(z-z_1)/(z_2-z_1);
> > dz_du = -tan(beta);
> >
> > function [gstop,isterminal,direction]=cross_z_z_min(u,z,data)
> > gstop = z-data.z_2;
> > isterminal = 1;
> > direction = [];
> >
> > 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?
> - - - - - - - - - -
> This differential equation can be solved without using 'ode15'. The solution (if I understand your code correctly) is:
>
> z = 1/(b2-b1)*(b2-asin(sin(b1)*exp((b2-b1)*u)))
>
> 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-(b2-b1)*z gets past pi/2.
>
> 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.
>
> Roger Stafford

----------------------------------------

thank you Roger for your reply.

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.
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.

I am not shure if i understood correctly your suggestion. should i split up the function, one part for
beta < pi /2
and one for
beta > pi/2?

Subject: ode, tan, singularity?

From: Roger Stafford

Date: 9 Mar, 2012 20:22:15

Message: 4 of 7

"Roland " <burgmann@gmx.de> wrote in message <jjch18$ens$1@newscl01ah.mathworks.com>...
> .........
> I am not shure if i understood correctly your suggestion. should i split up the function, one part for
> beta < pi /2
> and one for
> beta > pi/2?
- - - - - - - -
  When (b2-b1)/(z2-z1) < 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.

  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 (b2-b1)/(z2-z1) > 0 these trajectories are reversed in the u direction with asymptotes occurring as u approaches plus infinity.
  
Roger Stafford

Subject: ode, tan, singularity?

From: Roland

Date: 10 Mar, 2012 09:49:12

Message: 5 of 7

"Roger Stafford" wrote in message <jjdopn$ptc$1@newscl01ah.mathworks.com>...
> "Roland " <burgmann@gmx.de> wrote in message <jjch18$ens$1@newscl01ah.mathworks.com>...
> > .........
> > I am not shure if i understood correctly your suggestion. should i split up the function, one part for
> > beta < pi /2
> > and one for
> > beta > pi/2?
> - - - - - - - -
> When (b2-b1)/(z2-z1) < 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.
>
> 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 (b2-b1)/(z2-z1) > 0 these trajectories are reversed in the u direction with asymptotes occurring as u approaches plus infinity.
>
> Roger Stafford

------------------------------------

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.

Roland

Subject: ode, tan, singularity?

From: Roger Stafford

Date: 10 Mar, 2012 14:35:12

Message: 6 of 7

"Roland " <burgmann@gmx.de> wrote in message <jjf82o$g3m$1@newscl01ah.mathworks.com>...
> 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.
- - - - - - - - - -
  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.

  Probably it would be preferable to allow arclength, s, to be the independent variable and send two simultaneous equations to ode45: either

 beta = .... % linear or otherwise in z
 du/ds = cos(beta)
 dz/ds = -sin(beta)

or else

 beta = ....
 du/ds = -cos(beta)
 dz/ds = sin(beta)

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.)

Roger Stafford

Subject: ode, tan, singularity?

From: Roland

Date: 10 Mar, 2012 16:49:21

Message: 7 of 7

"Roger Stafford" wrote in message <jjfor0$3r6$1@newscl01ah.mathworks.com>...
> "Roland " <burgmann@gmx.de> wrote in message <jjf82o$g3m$1@newscl01ah.mathworks.com>...
> > 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.
> - - - - - - - - - -
> 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.
>
> Probably it would be preferable to allow arclength, s, to be the independent variable and send two simultaneous equations to ode45: either
>
> beta = .... % linear or otherwise in z
> du/ds = cos(beta)
> dz/ds = -sin(beta)
>
> or else
>
> beta = ....
> du/ds = -cos(beta)
> dz/ds = sin(beta)
>
> 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.)
>
> Roger Stafford

--------------------

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.
Thanks a lot for you help!

Tags for this Thread

What are tags?

A tag is like a keyword or category label associated with each thread. Tags make it easier for you to find threads of interest.

Anyone can tag a thread. Tags are public and visible to everyone.

Contact us