models/johnpye/zill.a4c

REQUIRE "ivpsystem.a4l";
REQUIRE "atoms.a4l";

(* 
	An example problem used in the text by Zill, 4th Ed, p 472.

	Dennis G Zill, A First Course in Differential Equations with Applications,
	PWS-KENTm Boston, 1989.

	Actually this is a single-variable integration and can be solved using the
	Runge-Kutta method (as illustrated by Zill).
*)
MODEL zill;
	x IS_A solver_var;
	y IS_A solver_var;
	dy_dx IS_A solver_var;

	xderiv: dy_dx = 2 * x * y;

METHODS

METHOD values;
	y := 1; x := 1; (* initial conditions *)
	dy_dx := 0;
END values;

METHOD specify;
	FIX x, y;
END specify;

METHOD ode_init;
	FREE x, y, dy_dx;

	x.ode_type := -1;
	dy_dx.ode_id := 1;   y.ode_id := 1; 
	dy_dx.ode_type := 2; y.ode_type := 1;

	x.obs_id := 1;
	y.obs_id := 2;
END ode_init;
	
METHOD on_load;
	RUN reset; RUN values;
	RUN ode_init;
END on_load;

METHOD self_test;
	(* assumes integration up to x = 1.5 *)
	ASSERT abs(x - 1.5) < 1e-7;
	ASSERT abs(y - 3.4904) < 0.001; (* this would be true if we used RK4, with h = 0.1 *)
END self_test;

END zill;
1	REQUIRE "ivpsystem.a4l";
2	REQUIRE "atoms.a4l";
3
4	(*
5	An example problem used in the text by Zill, 4th Ed, p 472.
6
7	Dennis G Zill, A First Course in Differential Equations with Applications,
8	PWS-KENTm Boston, 1989.
9
10	Actually this is a single-variable integration and can be solved using the
11	Runge-Kutta method (as illustrated by Zill).
12	*)
13	MODEL zill;
14	x IS_A solver_var;
15	y IS_A solver_var;
16	dy_dx IS_A solver_var;
17
18	xderiv: dy_dx = 2 * x * y;
19
20	METHODS
21
22	METHOD values;
23	y := 1; x := 1; (* initial conditions *)
24	dy_dx := 0;
25	END values;
26
27	METHOD specify;
28	FIX x, y;
29	END specify;
30
31	METHOD ode_init;
32	FREE x, y, dy_dx;
33
34	x.ode_type := -1;
35	dy_dx.ode_id := 1; y.ode_id := 1;
36	dy_dx.ode_type := 2; y.ode_type := 1;
37
38	x.obs_id := 1;
39	y.obs_id := 2;
40	END ode_init;
41
42	METHOD on_load;
43	RUN reset; RUN values;
44	RUN ode_init;
45	END on_load;
46
47	METHOD self_test;
48	(* assumes integration up to x = 1.5 *)
49	ASSERT abs(x - 1.5) < 1e-7;
50	ASSERT abs(y - 3.4904) < 0.001; (* this would be true if we used RK4, with h = 0.1 *)
51	END self_test;
52
53	END zill;