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 *)
END values;

METHOD specify;
	FIX x, y;
END specify;

METHOD ode_init;
	FREE x, y;
	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.0003; (* 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	END values;
25
26	METHOD specify;
27	FIX x, y;
28	END specify;
29
30	METHOD ode_init;
31	FREE x, y;
32	x.ode_type := -1;
33	dy_dx.ode_id := 1; y.ode_id := 1;
34	dy_dx.ode_type := 2; y.ode_type := 1;
35
36	x.obs_id := 1;
37	y.obs_id := 2;
38	END ode_init;
39
40	METHOD on_load;
41	RUN reset; RUN values;
42	RUN ode_init;
43	END on_load;
44
45	METHOD self_test;
46	(* assumes integration up to x = 1.5 *)
47	ASSERT abs(x - 1.5) < 1e-7;
48	ASSERT abs(y - 3.4904) < 0.0003; (* this would be true if we used RK4, with h = 0.1 *)
49	END self_test;
50
51	END zill;