runge_kutta.anubis
11.7 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
The Anubis/Paradize project.
Drawing solutions of first order differential equations
using the fourth order Runge-Kutta method.
Copyright (c) Alain Prouté 2005.
Author: Alain Prouté
The program defined in this file opens a host window. When a point of the window is
clicked upon, the solution of a given first order differential equation which passes
through this point is drawn.
First order differential equations (when solved with respect to y') have the following
form:
y' = f(x,y)
where y is the unknown function, x the free (time) variable, and f a function defined
on a subset of R^2. A solution of this equation is a function g (defined on some
interval of the reals), such that for all x in this interval, we have:
g'(x) = f(x,g(x))
The Cauchy-Lipschitz theorem asserts that if f is locally lipschitzian relative to its
second variable, then, by every point of the domain of f, passes one and exactly one
maximal (i.e. non extendable) solution of the equation.
Below is the function f used in this example (you may customize).
read tools/maybefloat.anubis
define Maybe(Float) _1 = success(1).
define Maybe(Float) _2 = success(2).
define Maybe(Float) _3 = success(3).
define Maybe(Float)
f
(
Maybe(Float) x,
Maybe(Float) y
) =
x - y*y.
abs(x)^y.
x+y*y*y.
x*x + y.
x*x - y*y.
x - y*x*y.
_3 * ((y*y)^(_1/_3)).
sin(x*y).
success(1)/y.
sin(y).
The Runge-Kutta method allows the computation of points on the graph of a solution
starting at some point of this graph. Starting at a point (x_0,y_0) in the domain of
f, we compute the points (x_n,y_n) belonging to the graph of the solution passing
through (x_0,y_0) as follows:
x_(n+1) = x_n + h
k_1 = hf(x_n,y_n)
k_2 = hf(x_n + h/2, y_n + k_1/2)
k_3 = hf(x_n + h/2, y_n + k_2/2)
k_4 = hf(x_n + h, y_n + k_3)
y_(n+1) = y_n + k_1/6 + k_2/3 + k_3/3 + k_4/6
The error made by this method against the actual solution is of the order of h^5 at
each step.
The next function computes y_(n+1) from x_n, y_n and h, using f.
define Maybe(Float)
runge_kutta_4
(
Maybe(Float) x_n,
Maybe(Float) y_n,
Maybe(Float) h,
(Maybe(Float) x, Maybe(Float) y) -> Maybe(Float) f
) =
with h2 = h/success(2),
k_1 = h * f(x_n,y_n),
k_2 = h * f(x_n + h2, y_n + k_1/success(2)),
k_3 = h * f(x_n + h2, y_n + k_2/success(2)),
k_4 = h * f(x_n + h, y_n + k_3),
y_n + ((k_1+k_4)/success(6)) + ((k_2+k_3)/success(3)).
A point of coordinates (x,y) must be painted into the buffer (of type HostImage). We
have two scaling numbers 'h_scale' and 'v_scale'.
define RGB white = rgb(255,255,255).
define Bool
paint
(
HostImage buffer,
Maybe(Float) x,
Maybe(Float) y,
Maybe(Float) h_scale,
Maybe(Float) v_scale,
Int32 xmax,
Int32 ymax
) =
with x1 = x/h_scale, y1 = y/v_scale,
if x1 is success(x2)
then if y1 is success(y2)
then with x3 = integral_part(x2),
y3 = integral_part(y2),
if (0 =< x3 & x3 < xmax & 0 =< y3 & y3 < ymax)
then (paint_rectangle(buffer,rect(x3,y3,x3+1,y3+1),white); true)
else false
else false
else false.
Now, if we want to draw a solution we must iterate the above function.
define One
draw_solution
(
Maybe(Float) x,
Maybe(Float) y,
Maybe(Float) h,
(Maybe(Float) x, Maybe(Float) y) -> Maybe(Float) f,
Maybe(Float) h_scale,
Maybe(Float) v_scale,
Maybe(Float) x_trans,
Maybe(Float) y_trans,
HostImage buffer,
Int32 xmax,
Int32 ymax
) =
if y is
{
failure then unique, // overflow
success(_) then
if paint(buffer,
x,
success(integer_to_float(ymax))*v_scale - y,
h_scale,
v_scale,
xmax,
ymax)
then draw_solution(x+h,
runge_kutta_4(x,y,h,
(Maybe(Float) _x, Maybe(Float) _y) |->
f(x+x_trans,y+y_trans)),
h,
f,
h_scale,v_scale,
x_trans,y_trans,
buffer,
xmax,ymax)
else unique
}.
The next function makes the paint_method. It just copies the content of the buffer to
the window, but only within the given rectangle.
define (HostWindow hw, List(Rectangle) rs) -> One
make_paint_method
(
HostImage buffer,
Int32 x_trans,
Int32 y_trans,
Int32 xmax,
Int32 ymax,
Int32 hs,
Int32 vs
) =
with red = rgb(255,0,0),
(HostWindow hw, List(Rectangle) rs) |->
forget(map((Rectangle r) |-> paint_image(hw,r,0,0,buffer),rs)); // copy buffer to window within r
paint_rectangle(hw,rect(
0,
ymax+y_trans-1,
xmax,
ymax+y_trans
),red);
paint_rectangle(hw,rect(
-x_trans+hs,
ymax+y_trans-4,
-x_trans+hs+1,
ymax+y_trans+4
),red);
paint_rectangle(hw,rect(
-x_trans,
0,
-x_trans+1,
ymax
),red);
paint_rectangle(hw,rect(
-x_trans-4,
ymax+y_trans-vs-1,
-x_trans+4,
ymax+y_trans-vs
),red).
The event handler.
define One
handle_keys
(
KeyboardState ks,
KeyboardKey kk,
Var(Maybe(Float)) h_scale,
Var(Maybe(Float)) v_scale
) =
if kk is
{
unknown then unique,
char(c) then if c = '+' then (h_scale <- (*h_scale) * success(0.5);
v_scale <- (*v_scale) * success(0.5)) else
if c = '-' then (h_scale <- (*h_scale) * success(2);
v_scale <- (*v_scale) * success(2)) else
unique,
unicode(a,b) then unique,
special(s) then unique
}.
define (HostWindow hw, HostWindowEvent(One) e) -> List(Rectangle)
make_event_handler
(
Rectangle default,
HostImage buffer,
Var(Maybe(Float)) h_scale,
Var(Maybe(Float)) v_scale,
Int32 h,
Int32 x_trans,
Int32 y_trans,
Int32 xmax,
Int32 ymax
) =
(HostWindow hw, HostWindowEvent(One) e) |->
if e is
{
quit then [],
expose then [default]
pointer_entering then [],
pointer_leaving then [],
key_down(ks,kk) then handle_keys(ks,
kk,
h_scale,
v_scale); [default],
mouse_move(ks,x,y) then [],
mouse_click(ks,mc,x,y) then if mc is left_down
then (
draw_solution( // to the right (if h > 0)
success(integer_to_float(x))**h_scale,
success(integer_to_float(ymax-y))**v_scale,
success(integer_to_float(h))**h_scale,
f,
*h_scale,
*v_scale,
success(integer_to_float(x_trans))**h_scale,
success(integer_to_float(y_trans))**v_scale,
buffer,
xmax,ymax);
draw_solution( // to the left (if h > 0)
success(integer_to_float(x))**h_scale,
success(integer_to_float(ymax-y))**v_scale,
success(integer_to_float(-h))**h_scale,
f,
*h_scale,
*v_scale,
success(integer_to_float(x_trans))**h_scale,
success(integer_to_float(y_trans))**v_scale,
buffer,
xmax,ymax);
[default]
)
else [],
mouse_wheel(ks,d,x,y) then [],
tick then [],
repaint(l) then [],
specific(_) then []
}.
Finally, here is our program. It first creates the buffer, then opens the host
window. If the host window has been opened, it must be shown by 'show', otherwise it
remains invisible.
global define One
runge_kutta
(
List(String) args
) =
with iw = (Int32)800,
ih = (Int32)600,
default = rect(0,0,iw,ih),
hs = (Int32)60,
vs = (Int32)60,
step = (Int32)1,
x_trans = (Int32)-(iw>>1),
y_trans = (Int32)-(ih>>1),
with buffer = to_host_image(create_rgba_image(iw,ih,rgba(100,100,100,255)),1),
if open_host_window(rect(100,50,100+iw,50+ih),
"Runge Kutta Method",
managed,
make_paint_method(buffer,x_trans,y_trans,iw,ih,hs,vs),
make_event_handler(default,buffer,
var(success(1.0)/success(integer_to_float(hs))),
var(success(1.0)/success(integer_to_float(vs))),
step,
x_trans,y_trans,
iw,ih),
(List(HostWindowEvent(One)) l) |-> l) is
{
failure then print("Cannot open window.\n"),
success(hw) then show(hw)
}.