ODE 1D Demo: corrected display and lsodar warnings 41/14341/2
Pierre-Aime Agnel [Fri, 4 Apr 2014 13:20:31 +0000 (15:20 +0200)]
There was an error in the roots function of the ODE causing the calculation to go over y_max

Change-Id: I0bdf5da85a796b4f35e27b3d38e19322e3f82ee9

scilab/modules/differential_equations/demos/ode/ode_1dvect/ode_1dvect.dem.sce

index 1bf2477..c6df372 100644 (file)
@@ -5,7 +5,7 @@
 
 function demo_ode_1dvect()
 
-    oldln=lines();
+    oldln = lines();
     lines(0)
     text = ["Examples of ODE''s in 1 dimension"; ..
     "A trajectory is plotted by clicking on the"; ..
@@ -15,34 +15,36 @@ function demo_ode_1dvect()
     "You can start over by clicking on the LEFT button again"; ..
     "  or stop everything by clicking on the RIGHT button." ];
 
-    messagebox(text,"modal");
+    messagebox(text, "modal");
 
-    function yprim=f(t,y),yprim=y^2-t;endfunction
-    function z=g(t,y),z=[y-ymin;y+ymax];endfunction
+    function yprim = f(t, y), yprim = y^2-t, endfunction
+    function z = g(t, y), z = [y-ymin, ymax-y], endfunction
 
     tmin = -3;
     tmax =  5;
     ymin = -3;
-    ymax =  3;
+    ymax =  4;
 
     t = tmin:1:tmax;
     y = ymin:1:ymax;
 
     my_handle = scf(100001);
-    clf(my_handle,"reset");
+    clf(my_handle, "reset");
     demo_viewCode("ode_1dvect.dem.sce");
 
-    nt = size(t,"*");
-    ny = size(y,"*");
-    fx = ones(nt,ny);
-    fy = feval(t,y,f);
-    champ(t,y,fx,fy);
-    xlabel("t","fontsize",3)
-    ylabel("y","fontsize",3)
-    a=gca();
-    a.margins(3)=0.2;
+    nt = size(t, "*");
+    ny = size(y, "*");
+    fx = ones(nt, ny);
+    fy = feval(t, y, f);
+    champ(t, y, fx, fy);
+    xlabel("t", "fontsize", 3)
+    ylabel("y", "fontsize", 3)
+    a = gca();
+    a.margins(3) = 0.2;
+    a.tight_limits = "on"
+    a.data_bounds = [tmin, ymin; tmax, ymax]
     title([_("ODE 1D vector field")
-    "dy/dt=  y^2 - t"],"fontsize",3)
+    "dy/dt=  y^2 - t"], "fontsize", 3)
 
     oldt0 = 10*tmax;
     oldy0 = 10*ymax;
@@ -50,31 +52,31 @@ function demo_ode_1dvect()
     dy                 = 0.1;
 
     while(%T)
-        [b,t0,y0]=xclick();
-        if or(b==[2 5 -1000]) then break,end;
-        if or(b==[0 3])& tmin<t0 & t0<tmax & ymin<y0 & y0<ymax then
-            t1=t0:0.1:tmax;
-            [sol1,rd]=ode("root",y0,t0,t1,1.d-2,1.d-4,f,2,g);
-            xpoly(t1(1:size(sol1,"*"))',sol1');
-            p1=gce();p1.thickness=2;p1.foreground=5;
+        [b, t0, y0] = xclick();
+        if or(b == [2 5 -1000]) then break, end;
+        if or(b == [0 3])& tmin<t0 & t0<tmax & ymin<y0 & y0<ymax then
+            t1 = t0:0.1:tmax;
+            [sol1, rd] = ode("root", y0, t0, t1, 1.d-2, 1.d-4, f, 2, g);
+            xpoly(t1(1:size(sol1, "*"))', sol1');
+            p1 = gce(), p1.thickness = 2, p1.foreground = 5;
 
-            t2=t0:-0.1:tmin;
-            [sol2,rd]=ode("root",y0,t0,t2,1.d-2,1.d-4,f,2,g);
-            xpoly(t2(1:size(sol2,"*"))',sol2');
+            t2 = t0:-0.1:tmin;
+            [sol2, rd] = ode("root", y0, t0, t2, 1.d-2, 1.d-4, f, 2, g);
+            xpoly(t2(1:size(sol2, "*"))', sol2');
 
-            p2=gce();p2.thickness=2;p2.foreground=5;
-            rep=[t0,y0,-1];
-            while rep(3)==-1 then
-                rep=xgetmouse();
-                t0=rep(1); y0=rep(2);
-                if (tmin<t0 & t0<tmax & ymin<y0 & y0<ymax) & (abs(t0-oldt0)>=dt | abs(y0-oldy0)>=dy) then
-                    t1=t0:0.1:tmax;
-                    [sol1,rd]=ode("root",y0,t0,t1,1.d-2,1.d-4,f,2,g);
-                    t2=t0:-0.1:tmin;
-                    [sol2,rd]=ode("root",y0,t0,t2,1.d-2,1.d-4,f,2,g);
-                    p1.data=[t1(1:size(sol1,"*"))' sol1'];
-                    p2.data=[t2(1:size(sol2,"*"))' sol2'];
-                    oldt0=t0; oldy0=y0;
+            p2 = gce(), p2.thickness = 2, p2.foreground = 5;
+            rep = [t0, y0, -1];
+            while rep(3) == -1 then
+                rep = xgetmouse();
+                t0 = rep(1); y0 = rep(2);
+                if (tmin<t0 & t0<tmax & ymin<y0 & y0<ymax) & (abs(t0-oldt0)> = dt | abs(y0-oldy0)> = dy) then
+                    t1 = t0:0.1:tmax;
+                    [sol1, rd] = ode("root", y0, t0, t1, 1.d-2, 1.d-4, f, 2, g);
+                    t2 = t0:-0.1:tmin;
+                    [sol2, rd] = ode("root", y0, t0, t2, 1.d-2, 1.d-4, f, 2, g);
+                    p1.data = [t1(1:size(sol1, "*"))' sol1'];
+                    p2.data = [t2(1:size(sol2, "*"))' sol2'];
+                    oldt0 = t0; oldy0 = y0;
                 end
             end
         end