1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18 package org.hipparchus.samples.ode;
19
20 import org.hipparchus.ode.ODEIntegrator;
21 import org.hipparchus.ode.ODEState;
22 import org.hipparchus.ode.OrdinaryDifferentialEquation;
23 import org.hipparchus.ode.nonstiff.DormandPrince853Integrator;
24 import org.hipparchus.ode.sampling.StepNormalizer;
25
26 import java.io.File;
27 import java.io.IOException;
28 import java.io.PrintStream;
29 import java.nio.charset.StandardCharsets;
30 import java.util.Locale;
31
32
33
34 public class LorenzPlotter {
35
36
37 private double duration;
38
39
40 private double step;
41
42
43 private double sigma;
44
45
46 private double rho;
47
48
49 private double beta;
50
51
52 private File output;
53
54
55 private int width;
56
57
58 private int height;
59
60
61 private double viewXRot;
62
63
64 private double viewZRot;
65
66
67
68 private LorenzPlotter() {
69 duration = 125.0;
70 step = 0.002;
71 sigma = 10;
72 rho = 28;
73 beta = 8.0 / 3.0;
74 output = null;
75 width = 1000;
76 height = 1000;
77 viewXRot = 70;
78 viewZRot = 20;
79 }
80
81
82
83
84 public static void main(String[] args) {
85 final LorenzPlotter plotter = new LorenzPlotter();
86 try {
87 for (int i = 0; i < args.length; ++i) {
88 switch (args[i]) {
89 case "--help" :
90 usage(0);
91 break;
92 case "--duration" :
93 plotter.duration = Double.parseDouble(args[++i]);
94 break;
95 case "--step" :
96 plotter.step = Double.parseDouble(args[++i]);
97 break;
98 case "--sigma" :
99 plotter.sigma = Double.parseDouble(args[++i]);
100 break;
101 case "--rhon" :
102 plotter.rho = Double.parseDouble(args[++i]);
103 break;
104 case "--beta" :
105 plotter.beta = Double.parseDouble(args[++i]);
106 break;
107 case "--output-dir" :
108 plotter.output = new File(args[++i]);
109 if (!(plotter.output.exists() && plotter.output.isDirectory() && plotter.output.canWrite())) {
110 System.err.format(Locale.US, "cannot generate output file in %s%n",
111 plotter.output.getAbsolutePath());
112 System.exit(1);
113 }
114 break;
115 case "--width" :
116 plotter.width = Integer.parseInt(args[++i]);
117 break;
118 case "--height" :
119 plotter.height = Integer.parseInt(args[++i]);
120 break;
121 case "--view" :
122 plotter.viewXRot = Double.parseDouble(args[++i]);
123 plotter.viewZRot = Double.parseDouble(args[++i]);
124 break;
125 default :
126 usage(1);
127 }
128 }
129
130 plotter.plot();
131
132 } catch (IndexOutOfBoundsException iobe) {
133 usage(1);
134 } catch (IOException ioe) {
135 System.err.println(ioe.getLocalizedMessage());
136 System.exit(1);
137 }
138
139 System.exit(0);
140
141 }
142
143
144
145
146 private static void usage(final int status) {
147 System.err.println("usage: java org.hipparchus.samples.ode.LorenzPlotter" +
148 " [--help]" +
149 " [--duration duration] [--step step]" +
150 " [--sigma sigma] [--rho rho] [--beta beta]" +
151 " [--output-dir directory]" +
152 " [--view xRot zRot]");
153 System.exit(status);
154
155 }
156
157
158
159
160 public void plot() throws IOException {
161 final ProcessBuilder pb = new ProcessBuilder("gnuplot").
162 redirectOutput(ProcessBuilder.Redirect.INHERIT).
163 redirectError(ProcessBuilder.Redirect.INHERIT);
164 pb.environment().remove("XDG_SESSION_TYPE");
165 final Process gnuplot = pb.start();
166 try (PrintStream out = new PrintStream(gnuplot.getOutputStream(), false, StandardCharsets.UTF_8.name())) {
167 final File outputFile;
168 if (output == null) {
169 out.format(Locale.US, "set terminal qt size %d, %d title 'Lorenz plotter'%n", width, height);
170 outputFile = null;
171 } else {
172 out.format(Locale.US, "set terminal pngcairo size %d, %d%n", width, height);
173 outputFile = new File(output, "Lorenz-attractor.png");
174 out.format(Locale.US, "set output '%s'%n", outputFile.getAbsolutePath());
175 }
176 out.format(Locale.US, "set xlabel 'X'%n");
177 out.format(Locale.US, "set ylabel 'Y'%n");
178 out.format(Locale.US, "set zlabel 'Z'%n");
179 out.format(Locale.US, "set key off%n");
180 out.format(Locale.US, "unset colorbox%n");
181 out.format(Locale.US, "set palette model RGB defined (0 0.45 0.66 0.85,0.25 0.725 0.9 1, 0.33 0.85 0.95 1, 0.33 0.68 0.82 0.65,0.4 0.58 0.75 0.53,0.575 0.95 0.9 0.75, 0.8 0.66 0.52 0.32,0.85 0.66 0.6 0.5, 1 0.95 0.95 0.95)%n");
182 out.format(Locale.US, "$data <<EOD%n");
183
184 ODEIntegrator integrator = new DormandPrince853Integrator(1.0e-3, 10.0, 1.0e-12, 1.0e-12);
185 integrator.addStepHandler(new StepNormalizer(step, (state, isLast) ->
186 out.format(Locale.US, "%.6f %.3f %.3f %.3f%n",
187 state.getTime() / duration,
188 state.getCompleteState()[0],
189 state.getCompleteState()[1],
190 state.getCompleteState()[2])));
191 integrator.integrate(new LorenzOde(),
192 new ODEState(0, new double[] { -8, 8, rho - 1 }),
193 duration);
194 out.format(Locale.US, "EOD%n");
195 out.format(Locale.US, "set view %f, %f%n", viewXRot, viewZRot);
196 out.format(Locale.US, "splot $data using 2:3:4:1 with lines lc palette%n");
197 if (output == null) {
198 out.format(Locale.US, "pause mouse close%n");
199 } else {
200 System.out.format(Locale.US, "output written to %s%n",
201 outputFile.getAbsolutePath());
202 }
203 }
204 }
205
206
207 private class LorenzOde implements OrdinaryDifferentialEquation {
208
209
210 @Override
211 public int getDimension() {
212 return 3;
213 }
214
215
216 @Override
217 public double[] computeDerivatives(final double t, final double[] y) {
218 return new double[] {
219 sigma * (y[1] - y[0]),
220 y[0] * (rho - y[2]) - y[1],
221 y[0] * y[1] - beta * y[2]
222 };
223 }
224 }
225
226 }