| Home | Trees | Indices | Help |
|---|
|
|
This file contains functions useful for solving differential equations
which occur commonly in a 1st semester differential equations course.
* desolve -- Computes the "general solution" to a 1st or 2nd order
ODE via maxima.
* desolve_laplace -- Solves an ODE using laplace transforms via maxima.
Initials conditions are optional.
* desolve_system -- Solves any size system of 1st order odes using maxima.
Initials conditions are optional.
* eulers_method -- Approximate solution to a 1st order DE, presented as a table.
* eulers_method_2x2 -- Approximate solution to a 1st order system of DEs,
presented as a table.
* eulers_method_2x2_plot -- Plots the sequence of points obtained from
Euler's method.
AUTHORS: David Joyner (3-2006) -- Initial version of functions
Marshall Hampton (7-2007) -- Creation of Python module and testing
Other functions for solving DEs are given in functions/elementary.py.
|
|||
|
|||
|
|||
|
|||
|
|||
|
|||
|
|||
|
|||
line = type line? for help and examples
|
|||
maxima = Maxima
|
|||
|
|||
Solves a 1st or 2nd order linear ODE via maxima. Initials conditions
are not given.
INPUT:
de -- a lambda expression representing the ODE
(eg, de = "diff(f(x),x,2)=diff(f(x),x)+sin(x)")
vars -- a list of strings representing the variables
(eg, vars = ["x","y"], if x is the independent
variable and y is the dependent variable)
EXAMPLES:
sage: from sage.calculus.desolvers import desolve
sage: t = var('t')
sage: x = function('x', t)
sage: de = lambda y: diff(y,t) + y - 1
sage: desolve(de(x(t)),[x,t])
'%e^-t*(%e^t+%c)'
AUTHOR: David Joyner (1-2006)
|
Solves an ODE using laplace transforms via maxima. Initials conditions
are optional.
INPUT:
de -- a lambda expression representing the ODE
(eg, de = "diff(f(x),x,2)=diff(f(x),x)+sin(x)")
vars -- a list of strings representing the variables
(eg, vars = ["x","f"], if x is the independent
variable and f is the dependent variable)
ics -- a list of numbers representing initial conditions,
with symbols allowed which are represented by strings
(eg, f(0)=1, f'(0)=2 is ics = [0,1,2])
EXAMPLES:
sage: from sage.calculus.desolvers import desolve_laplace
sage: x = var('x')
sage: f = function('f', x)
sage: de = lambda y: diff(y,x,x) - 2*diff(y,x) + y
sage: desolve_laplace(de(f(x)),["x","f"])
"x*%e^x*(?%at('diff(f(x),x,1),x=0))-f(0)*x*%e^x+f(0)*%e^x"
sage: desolve_laplace(de(f(x)),["x","f"],[0,1,2])
'x*%e^x+%e^x'
WARNING:
The second SAGE command in the above example sets the values of f(0) and f'(0)
in Maxima, so subsequent ODEs involving these variables will have these initial conditions
automatically imposed.
AUTHOR: David Joyner (1-2006,8-2007)
|
Solves any size system of 1st order odes using maxima. Initials conditions
are optional.
INPUT:
de -- a list of strings representing the ODEs in maxima notation
(eg, de = "diff(f(x),x,2)=diff(f(x),x)+sin(x)")
vars -- a list of strings representing the variables
(eg, vars = ["t","x","y"], where t is the independent variable
and x,y the dependent variables)
ics -- a list of numbers representing initial conditions
(eg, x(0)=1, y(0)=2 is ics = [0,1,2])
WARNING:
The given ics sets the initial values of the dependent vars in maxima, so
subsequent ODEs involving these variables will have these initial conditions
automatically imposed.
EXAMPLES:
sage: from sage.calculus.desolvers import desolve_system
sage: t = var('t')
sage: x = function('x', t)
sage: y = function('y', t)
sage: de1 = lambda z: diff(z[0],t) + z[1] - 1
sage: de2 = lambda z: diff(z[1],t) - z[0] + 1
sage: des = [de1([x(t),y(t)]),de2([x(t),y(t)])]
sage: vars = ["t","x","y"]
sage: desolve_system(des,vars)
['(1-y(0))*sin(t)+(x(0)-1)*cos(t)+1', '(x(0)-1)*sin(t)+(y(0)-1)*cos(t)+1']
sage: ics = [0,1,-1]
sage: soln = desolve_system(des,vars,ics); soln
['2*sin(t)+1', '1-2*cos(t)']
sage: solnx = lambda s: RR(eval(soln[0].replace("t","s")))
sage: solnx(3)
1.28224001611973
sage: solny = lambda s: RR(eval(soln[1].replace("t","s")))
sage: P1 = plot([solnx,solny],0,1)
sage: P2 = parametric_plot((solnx,solny),0,1)
Now type show(P1), show(P2) to view these.
AUTHOR: David Joyner (3-2006, 8-2007)
|
This implements Euler's method for finding numerically the solution of the 1st order
ODE y' = f(x,y), y(a)=c. The "x" column of the table increments from
x0 to x1 by h (so (x1-x0)/h must be an integer). In the "y" column,
the new y-value equals the old y-value plus the corresponding entry in the
last column.
*For pedagogical purposes only.*
EXAMPLES:
sage: from sage.calculus.desolvers import eulers_method
sage: x,y = PolynomialRing(QQ,2,"xy").gens()
sage: eulers_method(5*x+y-5,0,1,1/2,1)
x y h*f(x,y)
0 1 -2
1/2 -1 -7/4
1 -11/4 -11/8
sage: x,y = PolynomialRing(QQ,2,"xy").gens()
sage: eulers_method(5*x+y-5,0,1,1/2,1,method="none")
[[0, 1], [1/2, -1], [1, -11/4], [3/2, -33/8]]
sage: RR = RealField(sci_not=0, prec=4, rnd='RNDU')
sage: x,y = PolynomialRing(RR,2,"xy").gens()
sage: eulers_method(5*x+y-5,0,1,1/2,1,method="None")
[[0, 1], [1/2, -1.0], [1, -2.7], [3/2, -4.0]]
sage: RR = RealField(sci_not=0, prec=4, rnd='RNDU')
sage: x,y=PolynomialRing(RR,2,"xy").gens()
sage: eulers_method(5*x+y-5,0,1,1/2,1)
x y h*f(x,y)
0 1 -2.0
1/2 -1.0 -1.7
1 -2.7 -1.3
sage: x,y=PolynomialRing(QQ,2,"xy").gens()
sage: eulers_method(5*x+y-5,1,1,1/3,2)
x y h*f(x,y)
1 1 1/3
4/3 4/3 1
5/3 7/3 17/9
2 38/9 83/27
sage: eulers_method(5*x+y-5,0,1,1/2,1,method="none")
[[0, 1], [1/2, -1], [1, -11/4], [3/2, -33/8]]
sage: pts = eulers_method(5*x+y-5,0,1,1/2,1,method="none")
sage: P1 = list_plot(pts)
sage: P2 = line(pts)
sage: (P1+P2).show()
AUTHOR: David Joyner
|
This implements Euler's method for finding numerically the solution of the 1st order
system of two ODEs
x' = f(t, x, y), x(t0)=x0.
y' = g(t, x, y), y(t0)=y0.
The "t" column of the table increments from t0 to t1 by
h (so (t1-t0)/h must be an integer). In the "x" column, the new x-value equals the
old x-value plus the corresponding entry in the next (third) column.
In the "y" column, the new y-value equals the old y-value plus the
corresponding entry in the next (last) column.
*For pedagogical purposes only.*
EXAMPLES:
sage: from sage.calculus.desolvers import eulers_method_2x2
sage: t, x, y = PolynomialRing(QQ,3,"txy").gens()
sage: f = x+y+t; g = x-y
sage: eulers_method_2x2(f,g, 0, 0, 0, 1/3, 1,method="none")
[[0, 0, 0], [1/3, 0, 0], [2/3, 1/9, 0], [1, 10/27, 1/27], [4/3, 68/81, 4/27]]
sage: eulers_method_2x2(f,g, 0, 0, 0, 1/3, 1)
t x h*f(t,x,y) y h*g(t,x,y)
0 0 0 0 0
1/3 0 1/9 0 0
2/3 1/9 7/27 0 1/27
1 10/27 38/81 1/27 1/9
sage: RR = RealField(sci_not=0, prec=4, rnd='RNDU')
sage: t,x,y=PolynomialRing(RR,3,"txy").gens()
sage: f = x+y+t; g = x-y
sage: eulers_method_2x2(f,g, 0, 0, 0, 1/3, 1)
t x h*f(t,x,y) y h*g(t,x,y)
0 0 0.00 0 0.00
1/3 0.00 0.13 0.00 0.00
2/3 0.13 0.29 0.00 0.043
1 0.41 0.57 0.043 0.15
To numerically approximate y(1), where (1+t^2)y''+y'-y=0, y(0)=1,y'(0)=-1,
using 4 steps of Euler's method, first convert to a system:
y1' = y2, y1(0)=1; y2' = (y1-y2)/(1+t^2), y2(0)=-1.
sage: RR = RealField(sci_not=0, prec=4, rnd='RNDU')
sage: t, x, y=PolynomialRing(RR,3,"txy").gens()
sage: f = y; g = (x-y)/(1+t^2)
sage: eulers_method_2x2(f,g, 0, 1, -1, 1/4, 1)
t x h*f(t,x,y) y h*g(t,x,y)
0 1 -0.25 -1 0.50
1/4 0.75 -0.12 -0.50 0.29
1/2 0.63 -0.054 -0.21 0.19
3/4 0.63 -0.0078 -0.031 0.11
1 0.63 0.020 0.079 0.071
To numerically approximate y(1), where y''+ty'+y=0, y(0)=1,y'(0)=0:
sage: t,x,y=PolynomialRing(RR,3,"txy").gens()
sage: f = y; g = -x-y*t
sage: eulers_method_2x2(f,g, 0, 1, 0, 1/4, 1)
t x h*f(t,x,y) y h*g(t,x,y)
0 1 0.00 0 -0.25
1/4 1.0 -0.062 -0.25 -0.23
1/2 0.94 -0.11 -0.46 -0.17
3/4 0.88 -0.15 -0.62 -0.10
1 0.75 -0.17 -0.68 -0.015
AUTHOR: David Joyner
|
This plots the soln in the rectangle
(xrange[0],xrange[1]) x (yrange[0],yrange[1])
and plots using Euler's method the numerical solution of
the 1st order ODEs x' = f(t,x,y), x(a)=x0, y' = g(t,x,y), y(a) = y0.
*For pedagogical purposes only.*
===
The following example plots the solution to theta''+sin(theta)=0,
theta(0)=3/4, theta'(0) = 0. Type P[0].show() to plot the solution,
(P[0]+P[1]).show() to plot (t,theta(t)) and (t,theta'(t)).
===
EXAMPLES:
sage: from sage.calculus.desolvers import eulers_method_2x2_plot
sage: f = lambda z : z[2]; g = lambda z : -sin(z[1])
sage: P = eulers_method_2x2_plot(f,g, 0.0, 0.75, 0.0, 0.1, 1.0)
|
| Home | Trees | Indices | Help |
|---|
| Generated by Epydoc 3.0beta1 on Thu Jul 17 04:23:24 2008 | http://epydoc.sourceforge.net |