Report about Homotopy Continuation

I solved f(x) = 0 with the homotopy h(t, x) = tf(x) + (1 − t)g(x) by integrating the differntial equation

\{ \begin{matrix}x′(t) & = −\frac{h_{t}(t, x)}{h_{x}(t, x)} \\ x(0) & = x_{0} \end{matrix}

using Euler’s method with stepsize Δt = 1/10. First, a Maple script computes the terms h_{t} and h_{x} appearing in the

differential equation. In our test run we choose

f(x) = cos x − x, g(x) = f(x) − f(x_{0}) and x_{0} = 1.

The script func.mpl follows:

 1 # func.mpl -- compute forcing term for the differential equation
 2 # used in the homotopy continuation method for solving f(x)=0
 3 restart;
 4 f:=x->cos(x)-x;

5

 6 x0:=1.0;
 7 getx0:=unapply(x0,x);
 8 g:=unapply(f(x)-f(x0),x);
 9 df:=unapply(diff(f(x),x),x);
10 dg:=unapply(diff(g(x),x),x);
11 h:=unapply(t*f(x)+(1-t)*g(x),[t,x]);
12 dhdx:=unapply(diff(h(t,x),x),[t,x]);
13 dhdt:=unapply(diff(h(t,x),t),[t,x]);
14 F:=unapply(-dhdt(t,x)/dhdx(t,x),[t,x]);
15 with(CodeGeneration,C);
16 C(F);
17 fname:="func.i";
18 fd:=open(fname,WRITE);
19 fprintf(fd,"/* func.c -- by Maple */\n");
20 close(fd);
21 myopt:=output=fname,
22  declare=[x::numeric,t::numeric];
23 C(F,myopt);
24 C(f,myopt);
25 C(g,myopt);
26 C(h,myopt);
27 C(df,myopt);
28 C(dg,myopt);
29 C(dhdt,myopt);
30 C(dhdx,myopt);
31 C(getx0,myopt);

The Maple script creates a file func.i containing the needed C subroutines for computing a solution to the differential equation.

 1 /* func.c -- by Maple */
 2 #include <math.h>

3

 4 double F(double t, double x)
 5 {
 6  return (0.4596976941e0 / (t * (-sin(x) - 0.1e1) + (0.1e1 -
 7  t) * (-sin(x) - 0.1e1)));
 8 }

9

10 #include <math.h>

11

12 double f(double x)

13 {

14  return (cos(x) - x);

15 }

16

17 #include <math.h>

18

19 double g(double x)

20 {

21  return (cos(x) - x + 0.4596976941e0);

22 }

23

24 #include <math.h>

25

26 double h(double t, double x)

27 {

28  return (t * (cos(x) - x) + (0.1e1 - t) * (cos(x) - x +
29  0.4596976941e0));

30 }

31

32 #include <math.h>

33

34 double df(double x)

35 {

36  return (-sin(x) - 0.1e1);

37 }

38

39 #include <math.h>

40

41 double dg(double x)

42 {

43  return (-sin(x) - 0.1e1);

44 }

45

46 double dhdt(double t, double x)

48  return (-0.4596976941e0);

49 }

50

51 #include <math.h>

52

53 double dhdx(double t, double x)

54 {

55  return (t * (-sin(x) - 0.1e1) + (0.1e1 - t) * (-sin(x) - 0.1e1));

56 }

57

58 double getx0(double x)

59 {

60  return (0.10e1);

61 }

The program topy.c includes the file func.i and then solves the differential equation using Euler’s method.

 1 /*
 2  topy.c -- Solve f(x)=0 by the homotopy method
 3  Written Septeber 29, 2012 by Eric Olson for Math 701
 4 */

5

 6 #include <stdio.h>
 7 #include "func.i"

8

 9 #define N 10

10

11 main(){
12  printf("# topy -- Solve f(x)=0 by the homotopy method\n"
13  "# Written September 29, 2012 for Math 701 by Eric Olson\n\n");

14

15  int i;
16  double dt=1.0/N;
17  double x=getx0(0.0);
18  printf("# %2s %8s %24s\n","n","t_n","x_n");
19  for(i=0;i<N;i++){
20  double t=i*dt;
21  printf("%4d %8.4f %24.15e\n",i,t,x);
22  x=x+dt*F(t,x);

23 }

24  printf("%4d %8.4f %24.15e\n",N,N*dt,x);
25  printf("\n# f(x(1))=%24.15e\n",f(x));
26  return 0;

27 }

The resulting output

# topy -- Solve f(x)=0 by the homotopy method
# Written September 29, 2012 for Math 701 by Eric Olson
#  n      t_n                      x_n
0   0.0000    1.000000000000000e+00
1   0.1000    9.750363867857545e-01
2   0.2000    9.498849897486028e-01
3   0.3000    9.245342274005630e-01
4   0.4000    8.989718975600242e-01
5   0.5000    8.731851168990508e-01
6   0.6000    8.471602536636458e-01
7   0.7000    8.208828525619060e-01
8   0.8000    7.943375506400665e-01
9   0.9000    7.675079827553092e-01
10   1.0000    7.403766749982386e-01
# f(x(1))=  -2.162156049025259e-03

shows that a solution to f(x) = 0 is x ≈ 0.7404. Output from topy.c is stored in the file x.dat and processed with Gnu-

plot to create the graph

x(t)  1 0.8 0.6 0.4 0.2 0 1 0.95 0.9 0.85 0.8 0.75 0.7

using the Gnuplot script

1 set terminal pslatex
2 set output "graph.tex"

3 set style data lines
4 plot "x.dat" using 2:3 lt 3 ti "$x(t)$"

Note that all text output in x.dat is preceded by a ‘#’ sign which Gnuplot treats as comments when processing the data.