Application Center - Maplesoft

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

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

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：

Drawing of vectors at each point 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 at each lattice point . 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 ?, a direction field by drawing vectors on each lattice points is made as follows. Drawing a direction field on isocline  Let us say we have drawn a curve . At any point on the curve,   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 , , , , and vectors on those isoclines are drawn below.  Adding isocline：： (If you put of the value in the box and press the button (add), the isocline with the value are added.)

Range for draw：,

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

Solution curve

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

Initial condition：, , , ,
(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.