Applications of Numerical Continuation 

ODE-IVP Approach 

Hakan Tiftikci
Turkey
hakan.tiftikci@yahoo.com.tr 

Introduction  

In this  application, numerical continuation is applied to surface-surface intersection problem and kinematic analysis of a common mechanism called fourbar. There are numerous approaches to solve Numerical Continuation problem [Eugene,Kurt]. Approach employed in this work is to convert (nonlinear) equation system to a differential system which in turn describes an ODE-IVP problem, where initial condition comes from any solution on the continuation determined.
 

Overview of Numerical Continuation 

Numerical continuation methods follow the solution of one-parameter-family multivariate equation system as the parameter is varied continuously (or updated by small amounts). Symbolically such system is described by In usual form, and there are equations in unknowns so that for a fixed values of parameter , system becomes square and determinate. As the parameter is varied, point(s)/solution(s) of the system of equations trace the locus of solution as a curve in . If then the system is either underdetermined or overdetermined and usual procedure may be applied in least-square sense (using Moore-Penrose pseudoinverse as in Newton iteration). 

 

Typical introductory illustration in literature involves solving a multivariate equation system in presence of another "simpler" equation system whose solution is available or easier to solve than the original. In that case, the so-called "homotopy" function defined by allows continuing to zero of  starting with available solution obtained from simpler function and varying parameter from 0 to 1. 

 

For detailed explanation of available methods see [Eugene,Kurt] 

Employed Method 

The method employed in this work converts continuation problem to an ODE-IVP problem. The method required that image space and domain space dimensions are equal (excluding the parameter ) so there are equations in unknowns , 

 

 

... 

 

Geometrically this set describes hypersurfaces in (that constrain the variables ) so that the locus is a curve in  

Total differential of the equation system is 

 

 

... 

 

denoting 

 

 

total differentials in vectorial notation become dot product of gradients and differential  

 

 

... 

 

Above conditions geometrically indicate that the tangent to solution must be orthogonal to each of gradients So if a direction (vector) can be found such that it is orthogonal to all gradients, then solution curve can be traced by following this direction. Procedure given below determines this orthogonal vector. 

 

Given independent vectors in dimensional space, there is a unique direction in dimensional space orthogonal to all vectors given by 

 

where denotes basis vectors used to decompose/project vectors but treated as a symbol in this determinant expression (general orthogonal vector concept may be related to multilinear algebra with wedge products and Hodge star producing better notation, but to make implementation clear here, matrix notation is used) For example, in 2D there is unique perpendicular direction to 1 vector given by 

 

and in 3D there is unique direction orthogonal to 2 vectors given by 

 

 

 

which is recognized to be cross product of vectors   

Thus direction of solution curve tangent  may be obtained from gradients by 

 

Since this direction vector is is parallel to tangent of solution curve  where is a suitable scaling scalar, continuation problem is equivalant to differential equation  

 

where is an arbitrary time-like parameter and need not to be interpretable, is scaling parameter and is vector orthogonal to all gradients Note that for some problems it is preferrable to have one of the variables as differentiation parameter so that one can observe other variables with respect to and this approach is commonly used in kinematic analysis of mechanisms . 

Initialization 

 

>

restart:

 

Procedure Definitions 

 

Basic Rotation Matrices 

Following functions define rotation matrices about X,Y,Z axis. This matrices will be used to orient surfaces by Euler angles 

>	Rx := proc(t) <<1,0,0>|<0,cos(t),sin(t)>|<0,-sin(t),cos(t)>>; end proc:

 

>	Ry := proc(t) <<cos(t),0,-sin(t)>|<0,1,0>|<sin(t),0,cos(t)>>; end proc:

 

>	Rz := proc(t) <<cos(t),sin(t),0>|<-sin(t),cos(t),0>|<0,0,1>>; end proc:

 

Application to Surface-Surface intersections 

In this section continuation technique described above is applied to determine intersection of two surfaces. To illustrate the application of the method, surfaces are represented both in explicit and implicit form and continuation is applied in explicit-explicit, implicit-implicit, implicit-explicit combinations of two surfaces. 

 

One of the surfaces if ellipsoid of radii in 3 axis given implicitly by  

 

and explicitly by parameters  

 

 

 

where are angle of position vector with z-axis and xy-plane 

 

Other surface is cylinder of radius given implicitly by 

 

and explicitly by parameters  

 

 

 

If both surface representation are selected as explicit then equations take the form 

 

 

 

which is 3 scalar equations in four variables For both implicit representation there are two equations 

 

 

 

involving three variables  

 

For implicit-explicit case there are two sub alternatives. Equations are either 

 

 

 

having variables either 

or  

 

Following table summarize the configurations to represent two surfaces whose intersection is to be determined. Note that #Equaions=#Variables-1 for all configurations as required. 

 

Configuration	Equations	Variables	#Equations	#Variables
Explicit-Explicit			3	4
Implicit-Implicit	{		2	3
Explicit-Implicit			4	5
Explicit-Implicit			4	5

 

 

Summary of the application of the technique in following sections is described in following steps  

 

• ) setup equations in variables

 

• ) determine gradients by Jacobian of functions to yield

 

• ) stack with a row (symbolical) basis vectors

 

• ) determine orthogonal direction by determinant of

 

• ) replace each variable in orthogonal direction vector by its dependent form

 

• ) define the equation that assigns rates of variables to orthogonal direction

 

• ) determine a starting point (initial condition for ODE ) by fixing one of the parameters and solving for the rest (Newton's method, etc. can be used in numerical implementation)

 

• ) check solution graphically or comparing the residual to 0

 

 

Not: Steps 2-4 may be replaced by more concise implementation using Minors of matrix  

Define Surfaces 

 

Sphere implicit and Explicit (Parametric) Representations 

>	fSphere := (r, a1, a2, a3) -> r[1]^2/a1^2+r[2]^2/a2^2+r[3]^2/a3^2-1;

 

(3.1.1.1)

 

>	sSphere := (phi, theta, a1, a2, a3) -> <(a1*sin(phi)*cos(theta), a2*sin(phi)*sin(theta), a3*cos(phi))>;

 

(3.1.1.2)

 

 

Cylinder implicit and Explicit (Parametric) Representations 

>	fCylinder := (r, radius) -> r[1]^2+r[2]^2-radius^2;

 

(3.1.2.1)

 

>	sCylinder := (psi, h, radius) -> <radius*cos(psi), radius*sin(psi), h>;

 

(3.1.2.2)

 

 

Setup example surfaces 

 

Set radii of ellipsoid in X,Y,Z axis 

>	RS1 := 3; RS2 := 2; RS3 := 5/2;

 

 

 

	(3.1.3.1)

 

Set rotation matrix to orient ellipsoid (45 degrees turn around x-axis and then 60 degrees turn about rotated z-axis) 

>	C := Rx(Pi/4) . Rz(Pi/3);

 

(3.1.3.2)

 

Map global Cartesian coordinates to local frame  

>	R := LinearAlgebra:-Transpose(C) . <x, y, z>;

 

(3.1.3.3)

 

Determine implicit representation using transformed local coordinates 

>	fS := fSphere(R, RS1, RS2, RS3);

 

(3.1.3.4)

 

>	fS := expand(fS);

 

(3.1.3.5)

 

Determine explicit representation in local coordinates and transform it to global frame 

>	sS := C . sSphere(phi, theta, RS1, RS2, RS3);

 

(3.1.3.6)

 

Set cylinder radius 

>	RC := 3/2;

 

(3.1.3.7)

 

Determine implicit representation of cylinder 

>	fC := fCylinder([x, y, z], RC);

 

(3.1.3.8)

 

Determine explicit representation of cylinder 

>	sC := sCylinder(psi, h, RC);

 

(3.1.3.9)

 

 

 

 

 

>	ALI := .35: # Set ambient light intensity for rendering two surfaces

 

>	plotS := plot3d(convert(sS, 'list'), phi = 0 .. Pi, theta = 0 .. 2*Pi, axes = boxed, color = red, scaling = constrained, style = patch, glossiness = .9, light = ([45, 45, 1, 1, 1]), ambientlight = ([ALI, ALI, ALI])):

 

>	plotC := plot3d(convert(sC, 'list'), psi = 0 .. 2*Pi, h = -3.5 .. 3.5, axes = boxed, color = blue, scaling = constrained, style = patch, glossiness = .9, light = ([45, 45, 1, 1, 1]), ambientlight = ([ALI, ALI, ALI])):

 

>	plots[display]({plotC, plotS}, scaling = constrained, title = "two example surfaces : ellipsoid and cylinder");

 

 

 

 

Implicit-Implicit Case 

Define implicit equations for intersection of two surfaces. There are two scalar equations in three variables  

>	IIsystem := [fS, fC];

 

(3.2.1)

 

 

Compute Jacobian (~ gradients augmented) 

>	J := VectorCalculus[Jacobian](IIsystem, [x, y, z]);

 

(3.2.2)

 

Add basis vector symbols to compute perpendicular direction 

>	K := Matrix([[Matrix([ex, ey, ez])], [J]]);

 

(3.2.3)

 

>	OrthogonalVector := LinearAlgebra[Determinant](K);

 

(3.2.4)

 

equivalent computation using minor of Jacobian matrix 

>	seq((-1)^(i+1)*(LinearAlgebra[Minor])(K, 1, i), i = 1 .. 3);

 

(3.2.5)

 

>	OrthogonalVector, RHS := LinearAlgebra[GenerateMatrix]([OrthogonalVector], [ex, ey, ez]);

 

(3.2.6)

 

>	StateDerivatives := subs([x = x(t), y = y(t), z = z(t)], OrthogonalVector);

 

(3.2.7)

 

>

 

>	StateDerivatives := LinearAlgebra:-Row(StateDerivatives, 1);

 

(3.2.8)

 

>	LinearAlgebra:-Dimensions(%);

 

(3.2.9)

 

Define ODE-IVP description of continuation problem 

>	System := map(diff, <x(t), y(t), z(t)>, t)-LinearAlgebra:-Transpose(StateDerivatives);

 

(3.2.10)

 

To determine initial condition fix one of the variables   

>	IIsystemFixed := subs(y = 0, IIsystem);

 

(3.2.11)

 

and solve for others  

>	solution := solve(IIsystemFixed, {z, x});

 

(3.2.12)

 

>	allsolution := allvalues(solution[1]);

 

(3.2.13)

 

define initial condition 

>	IC := eval([x, y, z], `union`(allsolution[1], {y = 0}));

 

(3.2.14)

 

Define ODE-IVP  

>	System0 := [op(convert(System, 'list')), x(0) = IC[1], y(0) = IC[2], z(0) = IC[3]];

 

(3.2.15)

 

Solve numerically 

>	IIsolution := dsolve(System0, {x(t), y(t), z(t)}, 'numeric', method = rkf45, abserr = 0.1e-8, relerr = 0.1e-8);

 

(3.2.16)

 

plot variables against the ODE independent variable 

>	plots[odeplot](IIsolution, [t, x(t)], t = 0 .. 8, numpoints = 256,axes=frame);

 

 

>	plots[odeplot](IIsolution, [t, y(t)], t = 0 .. 8, numpoints = 256,axes=frame);

 

 

>	plots[odeplot](IIsolution, [t, z(t)], t = 0 .. 8, numpoints = 256,axes=frame);

 

 

plot solutions in the space of variables 

>	plots[odeplot](IIsolution, [x(t), y(t)], t = 0 .. 8, numpoints = 256, axes = boxed);

 

 

>	plots[odeplot](IIsolution, [x(t), z(t)], t = 0 .. 8, numpoints = 256, axes = boxed);

 

 

>	plots[odeplot](IIsolution, [y(t), z(t)], t = 0 .. 8, numpoints = 256, axes = boxed);

 

 

>	intersectionPlot := plots[odeplot](IIsolution, [x(t), y(t), z(t)], t = 0 .. 8, numpoints = 256, thickness = 3, color = green):

 

>	plots[display]({plotS, plotC, intersectionPlot});

 

 

Define the residuals 

>	resS := unapply(fS,[x,y,z])(x(t),y(t),z(t));

 

(3.2.17)

 

>	resC := unapply(fC,[x,y,z])(x(t),y(t),z(t));

 

(3.2.18)

 

Plot residuals 

>	plots[odeplot](IIsolution, [t,resS], t = 0 .. 8, numpoints = 256, axes = boxed);

 

 

>	plots[odeplot](IIsolution, [t,resC], t = 0 .. 8, numpoints = 256, axes = boxed);

 

 

>

 

Explicit-Explicit Case 

 

>	EEsystem := sS-sC;

 

(3.3.1)

 

>	J := VectorCalculus[Jacobian](EEsystem, [phi, theta, psi, h]);

 

(3.3.2)

 

>	K := Matrix([[Matrix([ephi, etheta, epsi, eh])], [J]]);

 

(3.3.3)

 

>	OrthogonalVector := LinearAlgebra[Determinant](K);

 

(3.3.4)

 

>	OrthogonalVector, RHS := LinearAlgebra[GenerateMatrix]([OrthogonalVector], [ephi, etheta, epsi, eh]);

 

(3.3.5)

 

>	StateDerivatives := subs([phi = phi(t), theta = theta(t), psi = psi(t), h = h(t)], OrthogonalVector);

 

(3.3.6)

 

>

 

>	StateDerivatives := LinearAlgebra:-Row(StateDerivatives, 1);

 

(3.3.7)

 

>

 

>	System := map(diff, <phi(t), theta(t), psi(t), h(t)>, t)-LinearAlgebra:-Transpose(StateDerivatives);

 

(3.3.8)

 

>

 

>	EEsystemFixed := subs(psi = 0, EEsystem);

 

(3.3.9)

 

>	solution := solve(convert(EEsystemFixed, 'list'), {theta, phi, h});

 

(3.3.10)

 

>	allsolution := allvalues(solution);

 

(3.3.11)

 

>	IC := eval([phi, theta, psi, h], `union`(allsolution[1], {psi = 0}));

 

(3.3.12)

 

>	evalf(IC);

 

(3.3.13)

 

>	System0 := [op(convert(System, 'list')), phi(0) = IC[1], theta(0) = IC[2], psi(0) = IC[3], h(0) = IC[4]];

 

(3.3.14)

 

>	EEsolution := dsolve(System0, {phi(t), theta(t), psi(t), h(t)}, 'numeric', method = rkf45, abserr = 0.1e-8, relerr = 0.1e-8);

 

(3.3.15)

 

>	rhs(EEsolution(1)[2]);

 

(3.3.16)

 

>	plots[odeplot](EEsolution, [t, h(t)], t = 0 .. 4, numpoints = 256,axes=frame);

 

 

>	plots[odeplot](EEsolution, [t, phi(t)], t = 0 .. 4, numpoints = 256,axes=frame);

 

 

>	plots[odeplot](EEsolution, [h(t), 180*phi(t)/Pi], t = 0 .. 4, numpoints = 256, axes = boxed);

 

 

>	plots[odeplot](EEsolution, [t, theta(t)], t = 0 .. 4, numpoints = 256,axes=frame);

 

 

>	plots[odeplot](EEsolution, [t, psi(t)], t = 0 .. 4, numpoints = 256,axes=frame);

 

 

>	plots[odeplot](EEsolution, [h(t), psi(t)], t = 0 .. 4, numpoints = 512, axes = boxed);

 

 

>	plots[odeplot](EEsolution, [phi(t), theta(t)], t = 0 .. 4, numpoints = 512, axes = boxed);

 

 

>

 

>	intersectionPlot := plots[spacecurve](convert(subs(psi = ''rhs(EEsolution(t)[4])'', h = ''rhs(EEsolution(t)[2])'', sC), 'list'), t = 0 .. 4, numpoints = 245, color = green, thickness = 3);

 

(3.3.17)

 

>	plots[display]({plotS, plotC, intersectionPlot}, title = "intersection of ellipsoid and from explicit (parametric) representations");

 

 

Define the residuals 

>	resSC := unapply(EEsystem,[phi,theta,psi,h])(phi(t),theta(t),psi(t),h(t));

 

(3.3.18)

 

>	plots[odeplot](EEsolution, [t,resSC[1]], t = 0 .. 4, numpoints = 512, axes = boxed); plots[odeplot](EEsolution, [t,resSC[2]], t = 0 .. 4, numpoints = 512, axes = boxed); plots[odeplot](EEsolution, [t,resSC[3]], t = 0 .. 4, numpoints = 512, axes = boxed);

 

 

 

 

 

 

>

 

>

 

Explicit-Implicit Case 

 

>	EIsystem1 := [op(convert(sS-<x,y,z>,'list')),fC];

 

(3.4.1)

 

>

 

>	EIsystem2 := [op(convert(sC-<x,y,z>,'list')),fS];

 

(3.4.2)

 

>	J1 := VectorCalculus[Jacobian](Vector[column](EIsystem1), [phi, theta, x,y,z]);

 

(3.4.3)

 

>	J2 := VectorCalculus[Jacobian](Vector[column](EIsystem2), [psi, h, x,y,z]);

 

(3.4.4)

 

>	K1 := Matrix([[Matrix([ephi, etheta, ex, ey,ez])], [J1]]);

 

(3.4.5)

 

>	K2 := Matrix([[Matrix([epsi, eh, ex, ey,ez])], [J2]]);

 

(3.4.6)

 

>	OrthogonalVector1 := LinearAlgebra[Determinant](K1);

 

(3.4.7)

 

>	OrthogonalVector2 := LinearAlgebra[Determinant](K2);

 

(3.4.8)

 

>	OrthogonalVector1, RHS1 := LinearAlgebra[GenerateMatrix]([OrthogonalVector1], [ephi, etheta, ex, ey,ez]);

 

(3.4.9)

 

>	OrthogonalVector2, RHS2 := LinearAlgebra[GenerateMatrix]([OrthogonalVector2], [epsi, eh, ex, ey,ez]);

 

(3.4.10)

 

>	StateDerivatives1 := subs([phi = phi(t), theta = theta(t), x = x(t), y = y(t), z=z(t)], OrthogonalVector1);

 

(3.4.11)

 

>	StateDerivatives2 := subs([psi = psi(t), h = h(t), x = x(t), y = y(t), z=z(t)], OrthogonalVector2);

 

(3.4.12)

 

>

 

>	StateDerivatives1 := LinearAlgebra:-Row(StateDerivatives1, 1);

 

(3.4.13)

 

>	StateDerivatives2 := LinearAlgebra:-Row(StateDerivatives2, 1);

 

(3.4.14)

 

>	System1 := map(diff, <phi(t), theta(t), x(t), y(t), z(t)>, t)-LinearAlgebra:-Transpose(StateDerivatives1);

 

(3.4.15)

 

>	System2 := map(diff, <psi(t), h(t), x(t), y(t), z(t)>, t)-LinearAlgebra:-Transpose(StateDerivatives2);

 

(3.4.16)

 

>	EIsystem2;

 

(3.4.17)

 

>	EIsystem1Fixed := eval(EIsystem1,[phi = Pi/6]);

 

(3.4.18)

 

>	EIsystem2Fixed := eval(EIsystem2,psi =0);

 

(3.4.19)

 

>	solution1 := solve(convert(EIsystem1Fixed, 'list'), {theta,x,y,z});

 

(3.4.20)

 

>	solution2 := solve(convert(EIsystem2Fixed, 'list'), {h,x,y,z});

 

(3.4.21)

 

>	allsolution1 := allvalues(solution1):

 

>	allsolution2 := allvalues(solution2);

 

(3.4.22)

 

>	IC1 := eval([phi, theta, x, y, z], `union`(allsolution1[1], {phi = Pi/6})):

 

>	IC2 := eval([psi, h, x, y, z], `union`(allsolution2[1], {psi = 0})):

 

>	evalf(IC1);

 

(3.4.23)

 

>	evalf(IC2);

 

(3.4.24)

 

>	System10 := [op(convert(System1, 'list')), phi(0) = IC1[1], theta(0) = IC1[2], x(0) = IC1[3], y(0) = IC1[4], z(0)=IC1[5]]:

 

>	System20 := [op(convert(System2, 'list')), psi(0) = IC2[1], h(0) = IC2[2], x(0) = IC2[3], y(0) = IC2[4], z(0)=IC2[5]]:

 

>	EI1solution := dsolve(System10, [phi(t), theta(t), x(t), y(t),z(t)], 'numeric', method = rkf45, abserr = 0.1e-8, relerr = 0.1e-8);

 

(3.4.25)

 

>	EI2solution := dsolve(System20, [psi(t), h(t), x(t), y(t),z(t)], 'numeric', method = rkf45, abserr = 0.1e-8, relerr = 0.1e-8);

 

(3.4.26)

 

>	rhs(EI1solution(1)[2]); rhs(EI2solution(1)[2]);

 

 

	(3.4.27)

 

>	plots[odeplot](EI1solution, [t, phi(t)], t = 0 .. 4, numpoints = 256);

 

 

>	plots[odeplot](EI1solution, [t, theta(t)], t = 0 .. 4, numpoints = 256);

 

 

>	plots[odeplot](EI1solution, [t, x(t)], t = 0 .. 4, numpoints = 256);

 

 

>	plots[odeplot](EI1solution, [t, y(t)], t = 0 .. 4, numpoints = 256);

 

 

>	plots[odeplot](EI1solution, [t, z(t)], t = 0 .. 4, numpoints = 256);

 

 

>	plots[odeplot](EI1solution, [phi(t)*180/Pi, 180*theta(t)/Pi], t = 0 .. 4, numpoints = 256, axes = boxed);

 

 

>

 

>	plots[odeplot](EI2solution, [psi(t)*180/Pi, h(t)], t = 0 .. 12, numpoints = 256, axes = boxed);

 

 

>	EI2solution(3);

 

(3.4.28)

 

>	intersectionPlot1 := plots[spacecurve](convert(subs(phi = ''rhs(EI1solution(t)[2])'', theta = ''rhs(EI1solution(t)[3])'', sS), 'list'), t = 0 .. 4, numpoints = 245, color = green, thickness = 3);

 

(3.4.29)

 

>	intersectionPlot2 := plots[spacecurve](convert(subs(psi = ''rhs(EI2solution(t)[2])'', h = ''rhs(EI2solution(t)[3])'', sC), 'list'), t = 0 .. 14, numpoints = 245, color = green, thickness = 3);

 

(3.4.30)

 

>	plots[display]({plotS, plotC, intersectionPlot1}, title = "intersection of ellipsoid and from explicit sphere - implicit cylinder representations");

 

 

>	plots[display]({plotS, plotC, intersectionPlot2}, title = "intersection of ellipsoid and from explicit cylinder - implicit sphere representations");

 

 

Define the residuals 

>	res1 := unapply(EIsystem1,[phi,theta,x,y,z])(phi(t),theta(t),x(t),y(t),z(t));

 

(3.4.31)

 

>	res2 := unapply(EIsystem2,[psi,h,x,y,z])(psi(t),h(t),x(t),y(t),z(t));

 

(3.4.32)

 

>

plots[odeplot](EI1solution, [t,res1[1]], t = 0 .. 4, numpoints = 512, axes = boxed); plots[odeplot](EI1solution, [t,res1[2]], t = 0 .. 4, numpoints = 512, axes = boxed); plots[odeplot](EI1solution, [t,res1[3]], t = 0 .. 4, numpoints = 512, axes = boxed); plots[odeplot](EI1solution, [t,res1[4]], t = 0 .. 4, numpoints = 512, axes = boxed);

 

 

 

 

 

>

plots[odeplot](EI2solution, [t,res2[1]], t = 0 .. 4, numpoints = 512, axes = boxed); plots[odeplot](EI2solution, [t,res2[2]], t = 0 .. 4, numpoints = 512, axes = boxed); plots[odeplot](EI2solution, [t,res2[3]], t = 0 .. 4, numpoints = 512, axes = boxed); plots[odeplot](EI2solution, [t,res2[4]], t = 0 .. 4, numpoints = 512, axes = boxed);

 

 

 

 

 

 

 

>

 

>

 

 

Application to Fourbar Mechanism 

A fourbar mechanism consists of 4 rigid links connected to each other by revolute joints (wikipedia). One technique is to solve.the so-called Loop-Closure-Equation. In this technique each link is represented by its length and orientation in fixed frame thus displacement between two successive revolute joints is given by which vanishes when summed for all links . For fourbar Loop-Closure-Equation (LCE) is given by following 

>	LCE := a1*exp(I*theta1)+a2*exp(I*theta2)+a3*exp(I*theta3)-a4;

 

(4.1)

 

When x and y components are evaluated 

>	LCEx := evalc(Re(LCE)); LCEy := evalc(Im(LCE));

 

 

	(4.2)

 

it is  noted that there are 2 equations and 3(=2+1) variables so continuation method described above can be applied to this problem  

>	LCE := [LCEx,LCEy]:

 

>	sol := solve({LCEx,LCEy},[theta2,theta3]):

 

>	whattype(sol); nops(sol);

 

 

	(4.3)

 

>	allsol := allvalues(sol):

 

>	whattype(allsol); nops([allsol]);

 

 

	(4.4)

 

>	J := VectorCalculus:-Jacobian(LCE, [theta1,theta2,theta3]);

 

(4.5)

 

>	K := Matrix([[Matrix([e1,e2,e3])],[J]]);

 

(4.6)

 

>	OrthogonalVector := LinearAlgebra[Determinant](K);

 

(4.7)

 

>	OrthogonalVector, RHS := LinearAlgebra[GenerateMatrix]([OrthogonalVector], [e1,e2,e3]);

 

(4.8)

 

>	StateDerivatives := subs([theta1 = theta1(t), theta2 = theta2(t), theta3 = theta3(t)], OrthogonalVector);

 

(4.9)

 

>	LinearAlgebra:-Dimensions(%);

 

(4.10)

 

>	StateDerivatives := LinearAlgebra:-Row(StateDerivatives, 1);

 

(4.11)

 

>	LinearAlgebra:-Dimensions(%);

 

(4.12)

 

>	System := map(diff, <theta1(t), theta2(t), theta3(t)>, t)-LinearAlgebra:-Transpose(StateDerivatives);

 

(4.13)

 

>

LCE;

 

(4.14)

 

>	kinematicDimensions := [a1=1,a2=2,a3=3/2,a4=2];

 

(4.15)

 

>	LCE0 := subs(kinematicDimensions,LCE);

 

(4.16)

 

>	eval(subs(theta1=0,LCE0));

 

(4.17)

 

>	sol := solve( eval(LCE0,[theta1=0]), {theta2,theta3});

 

(4.18)

 

>	IC := eval([theta1,theta2,theta3], `union`(sol, {theta1=0}));

 

(4.19)

 

>	System0 := [op(convert(subs(kinematicDimensions,System), 'list')), theta1(0) = IC[1], theta2(0) = IC[2], theta3(0) = IC[3]];

 

(4.20)

 

>	FourbarSolution := dsolve(System0, {theta1(t), theta2(t), theta3(t)}, 'numeric', method = rkf45, abserr = 0.1e-8, relerr = 0.1e-8);

 

(4.21)

 

>	FourbarSolution(1);

 

(4.22)

 

>	plots[odeplot](FourbarSolution, [t, theta1(t)], t = 0 .. 7, numpoints = 256);

 

 

>	plots[odeplot](FourbarSolution, [t, theta2(t)], t = 0 .. 7, numpoints = 256);

 

 

>	plots[odeplot](FourbarSolution, [t, theta3(t)], t = 0 .. 7, numpoints = 256);

 

 

>	plots[odeplot](FourbarSolution, [theta2(t), theta3(t)], t = 0 .. 7, numpoints = 256);

 

 

>	eval([theta1(t),theta2(t),theta3(t)], FourbarSolution(1));

 

(4.23)

 

>

plotFourbar := proc(KD,T) local t123,t1,t2,t3,X,Y,XY,sol; sol := FourbarSolution(T); t123 := eval([theta1(t),theta2(t),theta3(t)], sol); t1 := t123[1]; t2 := t123[2]; t3 := t123[3]; X := [0, a1*cos(t1), a1*cos(t1)+a2*cos(t2),a4]; Y := [0, a1*sin(t1), a1*sin(t1)+a2*sin(t2),0]; X := subs(KD, X); Y := subs(KD, Y); XY := zip( (a,b) -> [a,b], X,Y ); plots[listplot]( XY, color=blue ); end proc:

 

>	#plotFourbar(kinematicDimensions,0.2);

 

>	plots[animate]( plotFourbar, [kinematicDimensions ,T] ,T=0..5, frames=128, scaling=constrained,axes=boxed);

 

 

Define the residuals 

>	res := unapply(LCE0,[theta1,theta2,theta3])(theta1(t),theta2(t),theta3(t));

 

(4.24)

 

>	plots[odeplot](FourbarSolution, [t,res[1]], t = 0 .. 4, numpoints = 512, axes = boxed); plots[odeplot](FourbarSolution, [t,res[2]], t = 0 .. 4, numpoints = 512, axes = boxed);

 

 

 

>

 

Conclusions 

In this worksheet Numerical Continuation Method is applied to surface-surface intersection problem and kinematic analysis of mechanisms. For all 3 (explicit-explicit, explicit-implicit, implicit-implicit) representations of surfaces and LCE formulation of fourbar mechanism , method is verified by graphical means.  One important advantage of  the method described seems to be possibility of handling singular cases in a simple manner. In this respect, method is similar to Adjoint-Jacobian approach for (differential) inverse kinematics of serial manipulators. 

 

In literature, some other applications of numerical continuation methods can be found. Examples of application of the method include determination of trim point of an aircraft (for one selected varying parameter like speed, pitch angle, flight path angle, thrust, &c), envelope of moving planar shapes (for checking collision/interference of objects) and  inverse kinematics of serial manipulators.

Usual Numerical Continuation techniques include a corrector (step) to compensate for deviations from solution curve caused by finite steps taken by prediction step. A candidate to embed such correction to continous ODE-IVP approach is to define total residual scalar 

 

 

 

from which a direction for the minimization of residual is obtained by the gradient of , viz. 

 

 

 

so that nominal direction is bended towards the minimizing direction by  

 

 

 

with possibly a heurustic selection of suitable pairs 

References 

 

Eugene L. Allgower, Kurt Georg. Introduction to Numerical Continuation Methods, SIAM, 1987 

 

Numerical Continuation, Wikipedia 

 

Legal Notice: The copyright for this application is owned by the author(s). Neither Maplesoft nor the author are responsible for any errors contained within and are not liable for any damages resulting from the use of this material. This application is intended for non-commercial, non-profit use only. Contact the author for permission if you wish to use this application in for-profit activities.