[Scilab-users] Find the points of intersection of two curves

classic Classic list List threaded Threaded
12 messages Options
Hermes Hermes
Reply | Threaded
Open this post in threaded view
|

[Scilab-users] Find the points of intersection of two curves

Hello,
How could I improve the accuracy of intersection points?

clear;
funcprot(0);
timer();
function Sys=F(x)
x0=x(1);x1=x(2);
Sys(1)=x1^2-x0*x1+x0^2-1;
Sys(2)=sin(4*x1^2)+sin(5*x0^2);
endfunction  
function D=Draghilev(x)
x0=x(1);x1=x(2);x2=x(3);x00=ics(1);x01=ics(2);
D=[8*x01^2*x1*cos(4*x1^2)-8*x00*x01*x1*cos(4*x1^2)+8*x00^2*x1*cos(4*x1^2)-8*x1*cos(4*x1^2)-2*sin(4*x01^2)*x1-2*sin(5*x00^2)*x1+x0*sin(4*x01^2)+x0*sin(5*x00^2);
-sin(4*x01^2)*x1-sin(5*x00^2)*x1+2*x0*sin(4*x01^2)-10*x0*cos(5*x0^2)*x01^2+10*x0*cos(5*x0^2)*x00*x01+2*x0*sin(5*x00^2)-10*x0*cos(5*x0^2)*x00^2+10*x0*cos(5*x0^2);
-8*x1^2*cos(4*x1^2)+16*x0*x1*cos(4*x1^2)-20*x0*cos(5*x0^2)*x1+10*x0^2*cos(5*x0^2)];
endfunction

ics=[0;1]
ics=[ics;1];
disp(ics)
function [dxdt]=odes(t,x)
[dxdt]=Draghilev(x);
endfunction
N=100;
smin=0.0;
smax=5;
h=0.0001;

step=smax/N;
    t=smin:step:smax;
    t0=0;
    atol=h;
    LL= ode("fix",ics,t0,t,atol,odes)
disp(timer(),"tiempo de ejecucion C");
a=get("current_axes")//get the handle of the newly created axes
a.axes_visible="on"; // makes the axes visible
a.font_size=3; //set the tics label font size
a.x_location="middle"; //set the x axis position
plot(t,LL(1,1:$),"-cya",t,LL(2,1:$),"-g",t,LL(3,1:$),"-r")


s=LL(3,:)-LL(1,:);
ym=s(1:$-1).*s(2:$);
z=find(ym <= h);

yi=interpln([t;LL(1,:)],t(z+1));
plot(t(z+1),yi,"b*")

//
s=LL(3,:)-LL(2,:);
ym=s(1:$-1).*s(2:$);
z=find(ym <= 0);
I1=LL(2,find(ym <= h));
plot(t(z+1),I1,"g*")

Gracias




--
Sent from: http://mailinglists.scilab.org/Scilab-users-Mailing-Lists-Archives-f2602246.html
_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users
Christophe Dang Ngoc Chan Christophe Dang Ngoc Chan
Reply | Threaded
Open this post in threaded view
|

Re: {EXT} Find the points of intersection of two curves

Hello Hermes,

If I understand well,, the core of your question :

> De : Hermes
> Envoyé : vendredi 12 janvier 2018 10:47
>
> How could I improve the accuracy of intersection points?

Lies here:

> s=LL(3,:)-LL(1,:);
> ym=s(1:$-1).*s(2:$);
> z=find(ym <= h);

So you have 2 curves which y-values are LL(1,:) and LL(3,:)
and you search where the sign of the difference changes.
Correct?

So you suppose that both curves are continuous.

The refinement of the solution is more a mathematical topic than a Scilab one
and I'm afraid I don't have the skill and time to analyze your functions.

There are some general methods you could use:
use a polynomial interpolation instead of a linear interpolation,
i.e. choose an interval around the z values found above,
perform a polynomial regression
- you have an example of such code here : https://fr.wikibooks.org/wiki/D%C3%A9couvrir_Scilab/Calcul_num%C3%A9rique#R%C3%A9gression_polynomiale -
and search the local root of the polynomial.

You can also refine the step of calculation around the z values.

--
Christophe Dang Ngoc Chan
Mechanical calculation engineer
This e-mail may contain confidential and/or privileged information. If you are not the intended recipient (or have received this e-mail in error), please notify the sender immediately and destroy this e-mail. Any unauthorized copying, disclosure or distribution of the material in this e-mail is strictly forbidden.
_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users
Rafael Guerra Rafael Guerra
Reply | Threaded
Open this post in threaded view
|

Re: {EXT} Find the points of intersection of two curves

Hi,

This problem is equivalent to finding n-zeros over a given interval.
Matlab function fnzeros uses the following reference, which uses spline interpolation:

Mørken K. and Reimers, M. [2007] An unconditionally convergent method for computing zeros of splines and polynomials, Math. Comp., 76, 845--865

Regards,
Rafael

-----Original Message-----
From: users [mailto:[hidden email]] On Behalf Of Dang Ngoc Chan, Christophe
Sent: Friday, January 12, 2018 11:03 AM
To: Users mailing list for Scilab <[hidden email]>
Subject: Re: [Scilab-users] {EXT} Find the points of intersection of two curves

Hello Hermes,

If I understand well,, the core of your question :

> De : Hermes
> Envoyé : vendredi 12 janvier 2018 10:47
>
> How could I improve the accuracy of intersection points?

Lies here:

> s=LL(3,:)-LL(1,:);
> ym=s(1:$-1).*s(2:$);
> z=find(ym <= h);

So you have 2 curves which y-values are LL(1,:) and LL(3,:) and you search where the sign of the difference changes.
Correct?

So you suppose that both curves are continuous.

The refinement of the solution is more a mathematical topic than a Scilab one and I'm afraid I don't have the skill and time to analyze your functions.

There are some general methods you could use:
use a polynomial interpolation instead of a linear interpolation, i.e. choose an interval around the z values found above, perform a polynomial regression
- you have an example of such code here : https://fr.wikibooks.org/wiki/D%C3%A9couvrir_Scilab/Calcul_num%C3%A9rique#R%C3%A9gression_polynomiale - and search the local root of the polynomial.

You can also refine the step of calculation around the z values.

--
Christophe Dang Ngoc Chan
Mechanical calculation engineer
This e-mail may contain confidential and/or privileged information. If you are not the intended recipient (or have received this e-mail in error), please notify the sender immediately and destroy this e-mail. Any unauthorized copying, disclosure or distribution of the material in this e-mail is strictly forbidden.
_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users
_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users
Rafael Guerra Rafael Guerra
Reply | Threaded
Open this post in threaded view
|

Re: Find the points of intersection of two curves

In reply to this post by Hermes

Hi Hermes,

 

Here below a simple solution using linear interpolation that works fine with your data but would need some further refinement for general arbitrary inputs.

 

// START OF CODE

clear;

clf();

 

function Sys=F(x)

    x0=x(1);x1=x(2);

    Sys(1)=x1^2-x0*x1+x0^2-1;

    Sys(2)=sin(4*x1^2)+sin(5*x0^2);

endfunction

 

function D=Draghilev(x)

    x0=x(1);x1=x(2);x2=x(3);x00=ics(1);x01=ics(2);

    D=[8*x01^2*x1*cos(4*x1^2)-8*x00*x01*x1*cos(4*x1^2)+8*x00^2*x1*cos(4*x1^2)-8*x1*cos(4*x1^2)-2*sin(4*x01^2)*x1-2*sin(5*x00^2)*x1+x0*sin(4*x01^2)+x0*sin(5*x00^2);

    -sin(4*x01^2)*x1-sin(5*x00^2)*x1+2*x0*sin(4*x01^2)-10*x0*cos(5*x0^2)*x01^2+10*x0*cos(5*x0^2)*x00*x01+2*x0*sin(5*x00^2)-10*x0*cos(5*x0^2)*x00^2+10*x0*cos(5*x0^2);

    -8*x1^2*cos(4*x1^2)+16*x0*x1*cos(4*x1^2)-20*x0*cos(5*x0^2)*x1+10*x0^2*cos(5*x0^2)];

endfunction

 

function [dxdt]=odes(t, x)

    [dxdt]=Draghilev(x);

endfunction

 

ics=[0;1;1]

N=100;

smin=0.0;

smax=5;

h=0.0001;

 

step=smax/N;

t=smin:step:smax;

t0=0;

atol=h;

LL= ode("fix",ics,t0,t,atol,odes)

 

plot(t,LL(1,1:$),"-cya",t,LL(2,1:$),"-g",t,LL(3,1:$),"-r")

a=gca();

 

s= LL(3,:) - LL(1,:);

ym= s(1:$-1).*s(2:$);

z=find(ym <= 0);

t0 = t(z) - s(z).*(t(z+1)-t(z))./(s(z+1)-s(z));

y0 = interpln([t;LL(1,:)],t0);

plot(t0,y0,"b*")

 

s= LL(3,:) - LL(2,:);

ym= s(1:$-1).*s(2:$);

z=find(ym <= 0);

t0 = t(z) - s(z).*(t(z+1)-t(z))./(s(z+1)-s(z));

y0 = interpln([t;LL(2,:)],t0);

plot(t0,y0,"black*")

// END OF CODE

 

Regards,

Rafael

 

 

 


_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users
Rafael Guerra Rafael Guerra
Reply | Threaded
Open this post in threaded view
|

Re: Find the points of intersection of two curves

Hi Hermes,

Find here below a simple solution using linear interpolation that works fine with your data but would need some further refinement for general arbitrary inputs.

// START OF CODE
clear;
clf();

function Sys=F(x)
    x0=x(1);x1=x(2);
    Sys(1)=x1^2-x0*x1+x0^2-1;
    Sys(2)=sin(4*x1^2)+sin(5*x0^2);
endfunction

function D=Draghilev(x)
    x0=x(1);x1=x(2);x2=x(3);x00=ics(1);x01=ics(2);
    D=[8*x01^2*x1*cos(4*x1^2)-8*x00*x01*x1*cos(4*x1^2)+8*x00^2*x1*cos(4*x1^2)-8*x1*cos(4*x1^2)-2*sin(4*x01^2)*x1-2*sin(5*x00^2)*x1+x0*sin(4*x01^2)+x0*sin(5*x00^2);
    -sin(4*x01^2)*x1-sin(5*x00^2)*x1+2*x0*sin(4*x01^2)-10*x0*cos(5*x0^2)*x01^2+10*x0*cos(5*x0^2)*x00*x01+2*x0*sin(5*x00^2)-10*x0*cos(5*x0^2)*x00^2+10*x0*cos(5*x0^2);
    -8*x1^2*cos(4*x1^2)+16*x0*x1*cos(4*x1^2)-20*x0*cos(5*x0^2)*x1+10*x0^2*cos(5*x0^2)];
endfunction

function [dxdt]=odes(t, x)
    [dxdt]=Draghilev(x);
endfunction

ics=[0;1;1]
N=100;
smin=0.0;
smax=5;
h=0.0001;

step=smax/N;
t=smin:step:smax;
t0=0;
atol=h;
LL= ode("fix",ics,t0,t,atol,odes)

plot(t,LL(1,1:$),"-cya",t,LL(2,1:$),"-g",t,LL(3,1:$),"-r")
a=gca();

s= LL(3,:) - LL(1,:);
ym= s(1:$-1).*s(2:$);
z=find(ym <= 0);
t0 = t(z) - s(z).*(t(z+1)-t(z))./(s(z+1)-s(z));
y0 = interpln([t;LL(1,:)],t0);
plot(t0,y0,"b*")

s= LL(3,:) - LL(2,:);
ym= s(1:$-1).*s(2:$);
z=find(ym <= 0);
t0 = t(z) - s(z).*(t(z+1)-t(z))./(s(z+1)-s(z));
y0 = interpln([t;LL(2,:)],t0);
plot(t0,y0,"black*")
// END OF CODE

Regards,
Rafael


_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users
Hermes Hermes
Reply | Threaded
Open this post in threaded view
|

Re: Find the points of intersection of two curves

In reply to this post by Hermes
Hello Rafael,
with your proposal the result is as expected.

<http://mailinglists.scilab.org/file/t497622/GraphicScilabC.jpg>

Thank you
I have delayed responding because I wanted to arrive to the next step by my
means. But without results:
I want to achieve the following graph from the data obtained in LL.

<http://mailinglists.scilab.org/file/t497622/GraphicScilabD.jpg>

I have played with several configurations of the following graphical
function:
plot3d (LL (1, :), LL (2, :), 0 * ones (1,101), alpha = 35, theta = 45), for
the circle parallel to the xy plane.
plot3d (LL (1, :), LL (2, :), LL (3, :), alpha = 35, theta = 45), for the
function in space.
then it would be to put the intersection points.



--
Sent from: http://mailinglists.scilab.org/Scilab-users-Mailing-Lists-Archives-f2602246.html
_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users
Hermes Hermes
Reply | Threaded
Open this post in threaded view
|

Re: Find the points of intersection of two curves

The correct graph for the script shown should be the following:

<http://mailinglists.scilab.org/file/t497622/GraphicScilabD.jpg>

  my apologies



--
Sent from: http://mailinglists.scilab.org/Scilab-users-Mailing-Lists-Archives-f2602246.html
_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users
Rafael Guerra Rafael Guerra
Reply | Threaded
Open this post in threaded view
|

Re: Find the points of intersection of two curves

Hi Hermes,

Check this solution out:

// START OF CODE
clear;
clf();

function Sys=F(x)
    x0=x(1);x1=x(2);
    Sys(1)=x1^2-x0*x1+x0^2-1;
    Sys(2)=sin(4*x1^2)+sin(5*x0^2);
endfunction

function D=Draghilev(x)
    x0=x(1);x1=x(2);x2=x(3);x00=ics(1);x01=ics(2);
    D=[8*x01^2*x1*cos(4*x1^2)-8*x00*x01*x1*cos(4*x1^2)+8*x00^2*x1*cos(4*x1^2)-8*x1*cos(4*x1^2)-2*sin(4*x01^2)*x1-2*sin(5*x00^2)*x1+x0*sin(4*x01^2)+x0*sin(5*x00^2);
    -sin(4*x01^2)*x1-sin(5*x00^2)*x1+2*x0*sin(4*x01^2)-10*x0*cos(5*x0^2)*x01^2+10*x0*cos(5*x0^2)*x00*x01+2*x0*sin(5*x00^2)-10*x0*cos(5*x0^2)*x00^2+10*x0*cos(5*x0^2);
    -8*x1^2*cos(4*x1^2)+16*x0*x1*cos(4*x1^2)-20*x0*cos(5*x0^2)*x1+10*x0^2*cos(5*x0^2)];
endfunction

function [dxdt]=odes(t,x)
    [dxdt]=Draghilev(x);
endfunction

ics=[0;1;1]
N=100;
smin=0.0;
smax=5;
h=0.0001;

step=smax/N;
t=smin:step:smax;
t0=0;
atol=h;
LL= ode("fix",ics,t0,t,atol,odes)

s= LL(3,:);
ym= s(1:$-1).*s(2:$);
z= find(ym <= 0);
t0 = t(z) - s(z).*(t(z+1)-t(z))./(s(z+1)-s(z));
x0 = interpln([t;LL(1,:)],t0);
y0 = interpln([t;LL(2,:)],t0);

param3d(LL(1,:), LL(2,:), 0 * LL(1,:), alpha=35, theta=45); // circle in the xy plane
e = gce();
e.foreground =17;
param3d(LL(1,:), LL(2,:), LL(3,:), alpha=35, theta=45) // parametric function in 3d-space
param3d(x0, y0, 0*x0, alpha=35, theta=45) // zero-crossings
e = gce();
e.line_mode = "off";
e.mark_mode = "on"
e.mark_style = 3;
e.mark_foreground =5 ;
// END OF CODE

Regards,
Rafael
_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users
Rafael Guerra Rafael Guerra
Reply | Threaded
Open this post in threaded view
|

Re: Find the points of intersection of two curves

In previous email the code formatting got messed up, fixed it here below:

// START OF CODE
clear;
clf();

function Sys=F(x)
    x0=x(1);x1=x(2);
    Sys(1)=x1^2-x0*x1+x0^2-1;
    Sys(2)=sin(4*x1^2)+sin(5*x0^2);
endfunction

function D=Draghilev(x)
    x0=x(1);x1=x(2);x2=x(3);x00=ics(1);x01=ics(2);
    D=[8*x01^2*x1*cos(4*x1^2)-8*x00*x01*x1*cos(4*x1^2)+8*x00^2*x1*cos(4*x1^2)-8*x1*cos(4*x1^2)-2*sin(4*x01^2)*x1-2*sin(5*x00^2)*x1+x0*sin(4*x01^2)+x0*sin(5*x00^2);
    -sin(4*x01^2)*x1-sin(5*x00^2)*x1+2*x0*sin(4*x01^2)-10*x0*cos(5*x0^2)*x01^2+10*x0*cos(5*x0^2)*x00*x01+2*x0*sin(5*x00^2)-10*x0*cos(5*x0^2)*x00^2+10*x0*cos(5*x0^2);
    -8*x1^2*cos(4*x1^2)+16*x0*x1*cos(4*x1^2)-20*x0*cos(5*x0^2)*x1+10*x0^2*cos(5*x0^2)];
endfunction

function [dxdt]=odes(t,x)
    [dxdt]=Draghilev(x);
endfunction

ics=[0;1;1]
N=100;
smin=0.0;
smax=5;
h=0.0001;

step=smax/N;
t=smin:step:smax;
t0=0;
atol=h;
LL= ode("fix",ics,t0,t,atol,odes)

s= LL(3,:);
ym= s(1:$-1).*s(2:$);
z= find(ym <= 0);
t0 = t(z) - s(z).*(t(z+1)-t(z))./(s(z+1)-s(z));
x0 = interpln([t;LL(1,:)],t0);
y0 = interpln([t;LL(2,:)],t0);

param3d(LL(1,:), LL(2,:), 0 * LL(1,:), alpha=35, theta=45); // circle in the xy plane
e = gce();
e.foreground =17;
param3d(LL(1,:), LL(2,:), LL(3,:), alpha=35, theta=45); // parametric function in 3d-space
param3d(x0, y0, 0*x0, alpha=35, theta=45); // zero-crossings
e = gce();
e.line_mode = "off";
e.mark_mode = "on"
e.mark_style = 3;
e.mark_foreground =5 ;
// END OF CODE

Regards,
Rafael
_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users
Rafael Guerra Rafael Guerra
Reply | Threaded
Open this post in threaded view
|

Re: Find the points of intersection of two curves

(still problems with the cut and pasting of code from SciNotes, removed some of the comments now)

// START OF CODE
clear;
clf();

function Sys=F(x)
    x0=x(1);x1=x(2);
    Sys(1)=x1^2-x0*x1+x0^2-1;
    Sys(2)=sin(4*x1^2)+sin(5*x0^2);
endfunction

function D=Draghilev(x)
    x0=x(1);x1=x(2);x2=x(3);x00=ics(1);x01=ics(2);
    D=[8*x01^2*x1*cos(4*x1^2)-8*x00*x01*x1*cos(4*x1^2)+8*x00^2*x1*cos(4*x1^2)-8*x1*cos(4*x1^2)-2*sin(4*x01^2)*x1-2*sin(5*x00^2)*x1+x0*sin(4*x01^2)+x0*sin(5*x00^2);
    -sin(4*x01^2)*x1-sin(5*x00^2)*x1+2*x0*sin(4*x01^2)-10*x0*cos(5*x0^2)*x01^2+10*x0*cos(5*x0^2)*x00*x01+2*x0*sin(5*x00^2)-10*x0*cos(5*x0^2)*x00^2+10*x0*cos(5*x0^2);
    -8*x1^2*cos(4*x1^2)+16*x0*x1*cos(4*x1^2)-20*x0*cos(5*x0^2)*x1+10*x0^2*cos(5*x0^2)];
endfunction

function [dxdt]=odes(t,x)
    [dxdt]=Draghilev(x);
endfunction

ics=[0;1;1]
N=100;
smin=0.0;
smax=5;
h=0.0001;

step=smax/N;
t=smin:step:smax;
t0=0;
atol=h;
LL= ode("fix",ics,t0,t,atol,odes)

s= LL(3,:);
ym= s(1:$-1).*s(2:$);
z= find(ym <= 0);
t0 = t(z) - s(z).*(t(z+1)-t(z))./(s(z+1)-s(z));
x0 = interpln([t;LL(1,:)],t0);
y0 = interpln([t;LL(2,:)],t0);

param3d(LL(1,:), LL(2,:), 0 * LL(1,:), alpha=35, theta=45);
e = gce();
e.foreground =17;
param3d(LL(1,:), LL(2,:), LL(3,:), alpha=35, theta=45);
param3d(x0, y0, 0*x0, alpha=35, theta=45);
e = gce();
e.line_mode = "off";
e.mark_mode = "on";
e.mark_style = 3;
e.mark_foreground =5;
// END OF CODE

Regards,
Rafael

_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users
Hermes Hermes
Reply | Threaded
Open this post in threaded view
|

Re: Find the points of intersection of two curves


Hello Rafael,

<http://mailinglists.scilab.org/file/t497622/AllGraphicScilab.jpg>

Thankful for your help.
I think that in Scilab's help there should be examples of this type of
graphics.
Is there a site where you can see a gallery of "pretentious and intrepid"
graphics?
It would be a great help to assimilate the graphic potentialities of Scilab
to those who immigrated from other systems or software.




--
Sent from: http://mailinglists.scilab.org/Scilab-users-Mailing-Lists-Archives-f2602246.html
_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users
Rafael Guerra Rafael Guerra
Reply | Threaded
Open this post in threaded view
|

Re: Find the points of intersection of two curves

Hi Hermes,

You may want to check out this page in French which is very comprehensive:
https://perso.univ-rennes1.fr/philippe.roux/scilab/graphiques/fiche_graphiques.html

Regards,
Rafael

_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users