Application Center - Maplesoft
  Submit your work
Application Center Applications Geometrical View of Differential Equation dy/dx=f(x, y) Preview

App Preview:

Geometrical View of Differential Equation dy/dx=f(x, y)

You can switch back to the summary page by clicking here.

Learn about Maple
Download Application

 

Image 

Geometrical View of Differential Equation dy/dx=f(x, y) 

 

Yasuyuki Nakamura
Graduate School of Information Science, Nagoya University
A4-2(780), Furo-cho, Chikusa-ku, Nagoya, 464-8601, Japan
nakamura@nagoya-u.jp
http://www.phys.cs.is.nagoya-u.ac.jp/~nakamura/ 

Differential equation 

Let us consider a differential equation：Embedded component 

Drawing of vectors x, `/`(`*`(dy), `*`(dx)) at each point x, y is called a direction field.  

Below, we show two ways how to draw direction field.  

 

Drawing a dirrection field at each lattice point 

We draw a direction field x, `/`(`*`(dy), `*`(dx)) at each lattice point x, y. As we calculate each direction field one by one, this way is fit for drawing a direction field by a computer.  

When we make lattice points by dividing a reagion for drawing into Embedded component?Embedded component, a direction field by drawing vectors on each lattice points is made as follows. 

Draw of direction field 

Drawing a direction field on isocline 

Let us say we have drawn a curve f(x, y) = C. At any point on the curve,  `/`(`*`(dy), `*`(dx)) = C is held, which means any vector at the point on the curve has same direction. The cuve is called isocline. We can as many as vectors on the isocline and the way is good for drawing direction field by hands.  

Isocline (dotted blue line) with Embedded component, Embedded component, Embedded component, Embedded component, Embedded component and vectors on those isoclines are drawn below.isocline 

Adding isocline：Embedded component：add
(If you put of the value C in the box and press the button (add), the isocline with the value C are added.) 

 

Embedded component 

Embedded component 

 

Range for draw：Embedded componentEmbedded component,  Embedded componentEmbedded component 

Replot(If you change the plot range, please click [Replot] button.) 

Solution curve 

If the value y at x = x[0] (y(x[0]), initial condition) is given, the ODE can be solved. Solution curves are drawn with following initial condition.  

Initial condition：Embedded component, Embedded component, Embedded component, Embedded component, Embedded componentDraw
(Please put at least one initial condition.)  

 

Function 

 

> restart:
with(DocumentTools): with(DEtools): with(plots): with(StringTools): with(plottools):
p := []:
dp := []:
#ode := 1+x^2-y:
SetProperty('t_ode', 'value', 1+x^2-y):
SetProperty('grid_x', 'value', 20):
SetProperty('grid_y', 'value', 20):
SetProperty('x_min', 'value', 0):
SetProperty('x_max', 'value', 5):
SetProperty('y_min', 'value', -10):
SetProperty('y_max', 'value', 20):
SetProperty('t_isc_1', 'value', -6):
SetProperty('t_isc_2', 'value', -3):
SetProperty('t_isc_3', 'value', 0):
SetProperty('t_isc_4', 'value', 3):
SetProperty('t_isc_5', 'value', 6):
SetProperty('t_isc_a', 'value', -10):
SetProperty('t_iv_1', 'value', -10):
SetProperty('t_iv_2', 'value', 0):
SetProperty('t_iv_3', 'value', 5):
SetProperty('t_iv_4', 'value', 10):
SetProperty('t_iv_5', 'value', 15):
#ode := parse(GetProperty('t_ode', 'value')):
 

> set_ode := proc()
 global ode, ode_2;
 ode := parse(GetProperty('t_ode', 'value'));
 ode_2 := parse(StringTools[SubstituteAll](convert(ode, string), "y", "y(x)"));
end proc:
 

> clear_param := proc()
 global p, dp, l_c;
 p := [];
 dp := [];
 l_c := [];
end proc:
 

> set_range := proc()
 global x1, x2, y1, y2;
 x1 := parse(GetProperty('x_min', 'value'));
 x2 := parse(GetProperty('x_max', 'value'));
 y1 := parse(GetProperty('y_min', 'value'));
 y2 := parse(GetProperty('y_max', 'value'));
end proc:
 

> grid := proc()
 global ode, ode_2, x1, x2, y1, y2;
 local gx, gy;

 gx := parse(GetProperty('grid_x', 'value'));
 gy := parse(GetProperty('grid_y', 'value'));

 set_range();
 set_ode();
#  ode := parse(GetProperty('t_ode', 'value'));
 SetProperty('Plot0', 'value', DEplot(diff(y(x), x)=ode_2, y(x), x=x1..x2, y=y1..y2, dirfield=[gx, gy]));
end proc:
 

> iv := proc()
 global x1;
 local iv1, iv2, iv3, iv4, iv5, l_iv;

 l_iv := [];
 iv1 := DocumentTools[GetProperty](t_iv_1, value);
 iv2 := DocumentTools[GetProperty](t_iv_2, value);
 iv3 := DocumentTools[GetProperty](t_iv_3, value);
 iv4 := DocumentTools[GetProperty](t_iv_4, value);
 iv5 := DocumentTools[GetProperty](t_iv_5, value);

 if(iv1 <> "") then l_iv := [op(l_iv), y(x1)=parse(iv1)] end if;
 if(iv2 <> "") then l_iv := [op(l_iv), y(0)=parse(iv2)] end if;
 if(iv3 <> "") then l_iv := [op(l_iv), y(0)=parse(iv3)] end if;
 if(iv4 <> "") then l_iv := [op(l_iv), y(0)=parse(iv4)] end if;
 if(iv5 <> "") then l_iv := [op(l_iv), y(0)=parse(iv5)] end if;

 l_iv;
end proc:
 

> new_iso := proc()
 global ode, ode_2, x1, x2, y1, y2, l_c, p, dp;
 local c1, c2, c3, c4, c5, f, l, i, j, k, xi, yi, g;

#  ode := parse(GetProperty('t_ode', 'value'));
 set_range();
 set_ode();
 clear_param();

 c1 := GetProperty(t_isc_1, value);
 c2 := GetProperty(t_isc_2, value);
 c3 := GetProperty(t_isc_3, value);
 c4 := GetProperty(t_isc_4, value);
 c5 := GetProperty(t_isc_5, value);


 l_c := [];
 if(c1 <> "") then l_c := [op(l_c), parse(c1)] end if;
 if(c2 <> "") then l_c := [op(l_c), parse(c2)] end if;
 if(c3 <> "") then l_c := [op(l_c), parse(c3)] end if;
 if(c4 <> "") then l_c := [op(l_c), parse(c4)] end if;
 if(c5 <> "") then l_c := [op(l_c), parse(c5)] end if;

 iso();
end proc:
 

> iso := proc()
 global ode, ode_2, x1, x2, y1, y2, l_c, p, dp;
 local c1, c2, c3, c4, c5, f, l, i, j, k, xi, yi, g;

 if(IsSubSequence("y", convert(ode, string))) then

   xi := [seq(evalf((x2-x1)/20*i+x1), i=0..20)];

   if(IsSubSequence("x", SubstituteAll(convert(ode, string), "exp", ""))) then
   
     for i from 1 to nops(l_c) do
       f := [solve(ode=l_c[i], y)];
       l := [];
       for j from 1 to nops(f) do
         l := [op(l), f[j]];
         g := unapply(f[j], x);
         for k from 1 to 21 do
#            if(type(g(xi[k]), realcons)) then
           if( abs(Im(evalf(g(xi[k]))))<10^(-8) ) then
             dp := [op(dp), [xi[k], Re(evalf(g(xi[k])))]];
           end if;
         end do;
#          dp := [op(dp), seq([xi[k], g(xi[k])], k=1..21)];
       end do;
       p := [op(p), plot(l, x=x1..x2, y=y1..y2, color=blue, linestyle=dot)];
     end do;

   else

     for i from 1 to nops(l_c) do
       f := [my_fsolve(ode=l_c[i], y, y1, y2)];
       for j from 1 to nops(f) do
         p := [op(p), line([x1, f[j]], [x2, f[j]], color=blue, linestyle=dot)];
         for k from 1 to 21 do
           dp := [op(dp), [xi[k], f[j]]];
         end do;
       end do;
     end do;

   end if;
 else

   yi := [seq(evalf((y2-y1)/20*i+y1), i=0..20)];
   for i from 1 to nops(l_c) do
     f := [my_fsolve(ode=l_c[i], x, x1, x2)];
     for j from 1 to nops(f) do
       p := [op(p), line([f[j], y1], [f[j], y2], color=blue, linestyle=dot)];
       for k from 1 to 21 do
         dp := [op(dp), [f[j], yi[k]]];
       end do;
     end do;
   end do;

 end if;

 p := [op(p), DEplot(diff(y(x), x)=ode_2, y(x), x=x1..x2, y=y1..y2, dirfield=dp)];


 SetProperty('Plot1', 'value', display(seq(p[i], i=1..nops(p))));
end proc:
 

> add_iso := proc()
 global ode, ode_2, x1, x2, y1, y2, l_c, p, dp;
 local l, j, g, c_a, xi, yi, f, k, i;
 
 c_a := parse(GetProperty(t_isc_a, value));
 l_c := [op(l_c), c_a];


 if(IsSubSequence("y", convert(ode, string))) then

   xi := [seq(evalf((x2-x1)/20*i+x1), i=0..20)];

#    if(IsSubSequence("x", convert(ode, string))) then
   if(IsSubSequence("x", SubstituteAll(convert(ode, string), "exp", ""))) then
     f := [solve(ode=c_a, y)];
     l := [];
     for j from 1 to nops(f) do
       l := [op(l), f[j]];
       g := unapply(f[j], x);
       for k from 1 to 21 do
#          if(type(g(xi[k]), realcons)) then
         if( abs(Im(evalf(g(xi[k]))))<10^(-8) ) then
           dp := [op(dp), [xi[k], Re(evalf(g(xi[k])))]];
         end if;
       end do;
     end do;    
     p := [op(p[1..-2]), plot(l, x=x1..x2, y=y1..y2, color=blue, linestyle=dot)];


   else

     p := [op(p[1..-2])];
     f := [my_fsolve(ode=c_a, y, y1, y2)];
     for j from 1 to nops(f) do
       p := [op(p), line([x1, f[j]], [x2, f[j]], color=blue, linestyle=dot)];
       for k from 1 to 21 do
         dp := [op(dp), [xi[k], f[j]]];
       end do;
     end do;

   end if;  
   
 else

   yi := [seq(evalf((y2-y1)/20*i+y1), i=0..20)];
   p := [op(p[1..-2])];
   f := [my_fsolve(ode=c_a, x, x1, x2)];
   for j from 1 to nops(f) do
     p := [op(p), line([f[j], y1], [f[j], y2], color=blue, linestyle=dot)];
     for k from 1 to 21 do
       dp := [op(dp), [f[j], yi[k]]];
     end do;
   end do;

 end if;

 p := [op(p), DEplot(diff(y(x), x)=ode_2, y(x), x=x1..x2, y=y1..y2, dirfield=dp)];
 SetProperty('Plot1', 'value', display(seq(p[i], i=1..nops(p))));
end proc:
 

> orbit := proc()
 global ode, ode_2, x1, x2, y1, y2;
 local p2;
 set_range();

 SetProperty('Plot0', 'value', DEplot(diff(y(x), x)=ode_2, y(x), x=x1..x2, y=y1..y2, iv(), linecolor=red, animatecurves=true));
 SetProperty(Plot0, play, true);

 p2 := [op(p), DEplot(diff(y(x), x)=ode_2, y(x), x=x1..x2, y=y1..y2, iv(), linecolor=red, arrows=none)];
 SetProperty('Plot1', 'value', display(p2));
 SetProperty(Plot1, play, true);
end proc:
 

> my_fsolve := proc(expr, var, var_1, var_2)
 local sols, sol, is_exist, i;

 if( type(lhs(expr), polynom) ) then
   sols := fsolve(expr, var, var_1..var_2);
 else
   sols := {};
   is_exist := true;

   while( is_exist) do
     sol := fsolve(expr, var, avoid={op(sols)}, var=var_1..var_2);
     if (type(sol, realcons)) then
       sols := {op(sols), var=sol}
     else
       is_exist := false;
     end if;
   end do;

   sols := op(map(rhs, sols));
 end if;
 
 return sols;
end proc:
 

> re_plot := proc()
 global p, dp;

 set_range();
 p := [];
 dp := [];
 iso();
 grid();
end proc:
 

> grid(): new_iso():
 

>
 

 

Legal Notice: The copyright for this application is owned by the author(s). Neither Maplesoft nor the author are responsible for any errors contained within and are not liable for any damages resulting from the use of this material. This application is intended for non-commercial, non-profit use only. Contact the author for permission if you wish to use this application in for-profit activities.