ASKSAGE: Sage Q&A Forum - Latest question feedhttp://ask.sagemath.org/questions/Q&A Forum for SageenCopyright Sage, 2010. Some rights reserved under creative commons license.Sun, 16 Feb 2020 13:19:50 -0600I'm running across some strange syntax error that hasn't occured before I threw the cos function into the differential equation. Below is the code and the error:http://ask.sagemath.org/question/49954/im-running-across-some-strange-syntax-error-that-hasnt-occured-before-i-threw-the-cos-function-into-the-differential-equation-below-is-the-code-and/Here is the code
x,t = var('x,t')
plot1 = plot_slope_field(5*x-.3*x^2-2(1+cos(pi*t)),(t,0,10),(x,0,20),headaxislength=3,headlength=3,axes_labels=['$Time$','$Fish Pop.$'],title="No Harvesting",fontsize=13)
x = function('x')(t)
soln = lambda c:desolve_rk4(diff(x,t)==5*x-.3*x^2-2(1+cos(pi*t),x,ics=[0,c],end_points=[0,10],ivar=t,step=0.1)
plot2 = lambda c, color:list_plot(soln(c),plotjoined=True,axes_labels=['$Time$','$Fish Pop.$'],thickness=2,color=color)
show(plot1+sum(plot2(c,color=tuple(colormaps.hsv(1-c/13)[:3])) for c in [1,3.25,5.5,7.75,10]),ymin=0,ymax=20)
File "<ipython-input-27-1dec6aa982e7>", line 6
plot2 = lambda c, color:list_plot(soln(c),plotjoined=True,axes_labels=['$Time$','$Fish Pop.$'],thickness=Integer(2),color=color)
^
SyntaxError: invalid syntax
*The carrot was pointing at the "2" of "plot2"
Please help.PolarizedIceSun, 16 Feb 2020 13:19:50 -0600http://ask.sagemath.org/question/49954/I'm having trouble plotting this population differential equation with five particular solutions all on the same graph. Please help.http://ask.sagemath.org/question/49928/im-having-trouble-plotting-this-population-differential-equation-with-five-particular-solutions-all-on-the-same-graph-please-help/Here is some code
x,t = var('x,t')
plot1 = plot_slope_field(5*x-.3*x^2-0*t,(t,0,10),(x,0,20),headaxislength=3,headlength=3,axes_labels=['$Time$','$Fish Pop.$'],title="No Harvesting",fontsize=13)
x = function('x')(t)
soln = lambda c:desolve_rk4(diff(x,t)==5*x-.3*x^2-0*t,x,ics=[0,c],end_points=[0,10],ivar=t,step=0.1)
plot2 = lambda c:list_plot(soln(c),plotjoined=True,axes_labels=['$Time$','$Fish Pop.$'],thickness=2)
show(plot1+plot2(12),ymin=0,ymax=20)
Referencing the last line of code, I wish to input .1, 3, 6, 9, and 12 all at the same time.PolarizedIceFri, 14 Feb 2020 15:16:43 -0600http://ask.sagemath.org/question/49928/How to solve Sturm-Liouville problems in SageMath?http://ask.sagemath.org/question/47249/how-to-solve-sturm-liouville-problems-in-sagemath/I've conducted DuckDuckGo searches for "Sturm-Liouville sagemath" and received no usable results, and as a result I am here to ask, how can one solve Sturm-Liouville problems with SageMath? I know how to do so numerically with Julia, GNU Octave and Scilab, using matrices, but is there any way to solve at least some of them analytically in SageMath (e.g. $\frac{d^2 y}{dx^2}-xy = \lambda y$ has solutions $y_n=C\cdot \textit{Ai}(x+\lambda_n)$, where $\lambda_n $ are the zeros of the Airy $\textit{Ai}(x)$ function, if $y(0)=y(\infty)=0$)? Failing this, is there a way to numerically solve them in SageMath (without resorting to its interfaces with the languages I just mentioned)?
If it helps, here's the code I use in GNU Octave to solve the aforementioned SLE:
# Standard values are for transformation x=(y-L)/(y+L), y in 0, inf
# *lin are for x=2/L*y-1 (linear transformation)
clear all
format long e;
N=5000;
perc=0.266
Nfrag=round(perc*N);
n=0:N;
# 75 at N=1k; 162 for N=3k; 225 for N=5k;
L=226;
t=pi+n'*pi/N;
x=cos(t);
y=L*(cot(t/2)).^2;
ysub=y(2:N);
ylintrans = L/2*(x+1);
# Chebyshev polys of the first kind
T = cos(t*n);
# Now for arrays that do not include the endpoints
xsub = x(2:N);
Tsub = T(2:N,:);
Usub = diag(1./sqrt(1-xsub.^2))*sin(acos(xsub)*n);
dTsub = Usub*diag(n);
dT = [-(n.^2).*(-1).^(n); dTsub ; (n).^2];
dTL = diag(2*L./((y+L).^2))*dT;
dTlin = 2/L*dT;
D1 = dTL/T;
D1lin = dTlin/T;
D2 = D1^2;
D2lin = D1lin^2;
H = D2-diag(y);
Hlin = D2lin-diag(ylintrans);
H = H(2:N,2:N);
Hlin = Hlin(2:N,2:N);
[Y, Lam] = eig(H);
[Ylin, Lamlin] = eig(Hlin);
Lam = diag(Lam);
Lamlin = diag(Lamlin);
[Lam, IX] = sort(Lam, 'ascend');
[Lamlin, IXlin] = sort(Lamlin, 'ascend');
Y = Y(:,IX);
Ylin = Ylin(:,IXlin);
aizerochk = abs(airy(0,Lam));
aizerochklin = abs(airy(0,Lamlin));
L
errrms = sqrt(aizerochk(1:Nfrag)'*aizerochk(1:Nfrag)/(Nfrag+1))
errrmslin = sqrt(aizerochklin(1:Nfrag)'*aizerochklin(1:Nfrag)/(Nfrag+1))
ain=quad("airynorm",0,100);
figure(1)
plot(ysub(1:3*N/5), Y(1:3*N/5,1));
figure(2)
plot(ysub(1:3*N/5), Y(1:3*N/5,2));
figure(3)
plot(ysub(1:3*N/5), Y(1:3*N/5,Nfrag));
. It uses a Chebyshev spectral method on the extrema grid, and then analyses the results.Fusion809Wed, 24 Jul 2019 00:17:12 -0500http://ask.sagemath.org/question/47249/Solving an ODE and simplifying the resulthttp://ask.sagemath.org/question/46517/solving-an-ode-and-simplifying-the-result/ I'm interested in solving the differential equation $$3 h' + 3 h^2 = c_1,$$ where $c_1$ is a positive real number.
var('t')
var('c1', latex_name=r'c_1')
h = function('h')(t)
eq = -3*h^2 + c1 - 3*diff(h, t)
eq_sol = desolve(eq, h, ivar=t, contrib_ode=True)
The above code works, but it's not solved explicitly for $h$, so
h_sol = solve(eq_sol, h)
h_sol = h_sol[0]
h_sol
This gives something like $$h\left(t\right) = \frac{\sqrt{3} \sqrt{c_{1}} {\left(e^{\left(\frac{2}{3} \, \sqrt{3} C \sqrt{c_{1}} + \frac{2}{3} \, \sqrt{3} \sqrt{c_{1}} t\right)} + 1\right)}}{3 \, {\left(e^{\left(\frac{2}{3} \, \sqrt{3} C \sqrt{c_{1}} + \frac{2}{3} \, \sqrt{3} \sqrt{c_{1}} t\right)} - 1\right)}},$$
in sage notation (non-LaTeX) it starts like
h(t) == 1/3*sqrt(3)*sqrt(c1)* ...
**Question 1:** Is there a way to allocate to the solution (i.e. `h_sol`) the RHS of the above? without the `h(t) == ` part.
I had to set by hand (it is ease, but it would be nice to automatize the allocation)
var('C') # the integration constant introduced above
h_sol = 1/3*sqrt(3)*sqrt(c1)* ...
Then, by simply looking at the solution it is clear that it can be simplified. I tried things like
h_sol = h_sol.canonicalize_radical()
h_sol = h_sol.collect_common_factors()
h_sol = h_sol.simplify_rectform(complexity_measure = None)
but none of them returns the expected result, which could be obtained from Mathematica's kernel
mathematica("DSolve[3*h'[t] + 3*h[t]^2 == C[1], h[t], t]//FullSimplify")
$$ \sqrt{\frac{c_1}{3}} \tanh\left( \sqrt{\frac{c_1}{3}} (t - 3 c_2) \right) $$
**Question 2:** How could the expression `h_sol` be manipulated to obtain the hyperbolic tangent?DoxTue, 14 May 2019 02:41:24 -0500http://ask.sagemath.org/question/46517/Subtitute functions - in a differential equation - Sagemanifoldhttp://ask.sagemath.org/question/44557/subtitute-functions-in-a-differential-equation-sagemanifold/Dear community,
I have a differential equation that depends on a function $\xi(t)$, but is a component of a tensor (calculated with `sagemanifold`)
M = Manifold(4, 'M', latex_name=r"\mathcal{M}")
U.<t,r,th,ph> = M.chart(r't r:(0,+oo) th:(0,pi):\theta ph:(0,2*pi):\phi')
xi = function('xi')(t)
f = function('f')(t)
h = function('h')(t)
g = function('g')(t)
Ric = M.tensor_field(0,2, 'Ric')
Ric[0,0] = 3/2*f*h - 3/4*h^2 + 3/2*f*xi + 3/4*xi^2 - 3/2*diff(h, t) - 3/2*diff(xi, t)
Ric.display()
I'd like to define the restriction to $\xi = 0$, and assign it to a new tensor
Ric0 = M.tensor_field(0,2, 'Ric0')
Ric0[0,0] = Ric[0,0].substitute_expression({xi:0, diff(xi, t):0})
Ric0.display()
but I get an `AttributeError` because the
AttributeError: 'ChartFunctionRing_with_category.element_class' object has no attribute 'substitute_expression'
I know that it works for functions
var('t')
xi = function('xi')(t)
f = function('f')(t)
h = function('h')(t)
ode = 3/2*f*h - 3/4*h^2 + 3/2*f*xi + 3/4*xi^2 - 3/2*diff(h, t) - 3/2*diff(xi, t)
ode0 = ode.substitute_expression({xi:0, diff(xi, t):0})
<h2>Question</h2>
Is there a way to substitute functions that are not `ππππ.ππ’ππππππ.ππ‘ππππππππ.π΄π‘ππππππππ` but `ChartFunctionRing_with_category.element_class`?DoxMon, 03 Dec 2018 11:33:37 -0600http://ask.sagemath.org/question/44557/I'm trying to substitute a value in a Laplace transformhttp://ask.sagemath.org/question/44466/im-trying-to-substitute-a-value-in-a-laplace-transform/ Hello, I'm trying to solve a DE by using Laplace transforms, but when I get the lapalce transfrom of my DE I get this part s*D[0](y)(0) and I don't know how can I substitute D[0], I tried D[0](y)(0) == -1) but I get an error that says that D is not defined, I can subs y(0) == 1 with no issues. So my question is how can I substitute the value for y'(0)=-1 in the Laplace transform.
thank you for your help !.
t = var('t')
s = var('s')
y = function('y')(t)
ED = diff(y,t,2)-2*diff(y,t)-3*y==4
ED_Lap = ED.laplace(t,s)
ED_Lap = solve(ED.laplace(t,s),laplace(y(t), t, s))
print(ED_Lap)
show(ED_Lap[0].rhs().subs(y(0)==1,D[0](y)(0)==-1).inverse_laplace(s,t)) # I don't know how to subs DbemolMon, 26 Nov 2018 23:09:05 -0600http://ask.sagemath.org/question/44466/plotting autonomous differential equationshttp://ask.sagemath.org/question/43748/plotting-autonomous-differential-equations/I'm trying to plot the slope field and some solutions to
$$\frac{dx}{dt}=x^2-4x$$
I was able to get the slope field to plot with:
x,y=var('x','t')
plot_slope_field(x^2-4*x,(t,0,5),(x,-4,8))
I can't figure out how to the solutions to plot. I can get a general solution:
t = var('t')
x = function('x')(t)
f=desolve(diff(x,t) == x^2-4*x,[x,t])
show(f)
$$\newcommand{\Bold}[1]{\mathbf{#1}}\frac{1}{4} \, \log\left(x\left(t\right) - 4\right) - \frac{1}{4} \, \log\left(x\left(t\right)\right) = C + t$$
If I try it with initial values, say ics=[2,2], I get imaginary solutions:
$$\newcommand{\Bold}[1]{\mathbf{#1}}\frac{1}{4} \, \log\left(x\left(t\right) - 4\right) - \frac{1}{4} \, \log\left(x\left(t\right)\right) = \frac{1}{4} i \, \pi + t - 2$$
I was able to plot solutions with another windows tool, winplot, but really want to do it in sage. Any clues on how I can do that?AndyHTue, 25 Sep 2018 02:00:20 -0500http://ask.sagemath.org/question/43748/Solving a differential equation - From SageManifoldhttp://ask.sagemath.org/question/43370/solving-a-differential-equation-from-sagemanifold/ I'm interested in solving a differential equation obtained from a calculation from `sagemanifold`. I'll give a simple example, but the reason is that in general the equations **are not** as simple and one can introduce errors copying the equations.
I'd like to solve the equation given by the Ricci flat condition (for a certain affine connection).
reset()
M = Manifold(4, 'M', latex_name=r"\mathcal{M}")
U.<t,r,th,ph> = M.chart(r't r:(0,+oo) th:(0,pi):\theta ph:(0,2*pi):\phi')
nab = M.affine_connection('nabla', r'\nabla'); nab
k = var('k', latex_name=r'\kappa')
s = sqrt(1 - k * r**2)
f = function('f')(t)
g = function('g')(t)
h = function('h')(t)
nab[0,0,0] = f
nab[0,1,1] = g / (1 - k * r**2)
nab[0,2,2] = r**2 * g
nab[0,3,3] = r**2 * sin(th)**2 * g
nab[1,0,1] = h
nab[1,1,0] = h
nab[1,1,1] = k * r / (1 - k * r**2)
nab[1,2,2] = k * r**3 - r
nab[1,3,3] = (k * r **3 - r) * sin(th)**2
nab[2,0,2] = h
nab[2,1,2] = 1 / r
nab[2,2,0] = h
nab[2,2,1] = 1 / r
nab[2,3,3] = - cos(th) * sin(th)
nab[3,0,3] = h
nab[3,1,3] = 1 / r
nab[3,2,3] = cos(th) / sin(th)
nab[3,3,0] = h
nab[3,3,1] = 1 / r
nab[3,3,2] = cos(th) / sin(th)
nab.display()
Ric = M.tensor_field(0,2, 'R', latex_name=r'R')
Ric = nab.ricci()
print("Ricci tensor")
Ric.display_comp()
Now, the `Ric[0,0]` component yield a simple equation, and can be solved by the command
desolve(diff(h,t) + h**2 - f*h, h, ivar=t, contrib_ode=True)
However, I'd like to be able of telling *Sage* the following
desolve( Ric[0,0], h, ivar=t, contrib_ode=True)
but it seems that there is an incompatibility of types here...
type(diff(h,t) + h**2 - f*h)
returns `<ππ’ππ'ππππ.ππ’ππππππ.ππ‘ππππππππ.π΄π‘ππππππππ'>`, while
type(Ric[0,0])
returns `<πππππ'ππππ.πππππππππ.πππππβ―ππππ.π²πππππ΅πππππππππππβ―π πππβ―ππππππππ’.πππππππβ―πππππ'>`
------
**Question:**
What can I do to use the results from `sagemanifolds` to solve the differential equation?DoxMon, 13 Aug 2018 13:05:22 -0500http://ask.sagemath.org/question/43370/How to solve and plot y' = x^2 + y^3http://ask.sagemath.org/question/43369/how-to-solve-and-plot-y-x2-y3/Hi.
My first post!
I am trying to complete this exercise:
12. Although it might not be obvious from the differential equation, its solution could βbehave badlyβ near a point x at which we wish to approximate y(x). Numerical procedures may give widely differing results near this point. Let y(x) be the solution of the initial-value problem y' = x^2 + y^3, y(1) = 1.
(a) Use a numerical solver to graph the solution on the interval [1, 1.4].
(b) Using the step size h = 0.1, compare the results obtained from Eulerβs method with the results from the improved Eulerβs method in the approximation of y(1.4).
Please help?
THANK YOU!RobertWebbMon, 13 Aug 2018 12:41:49 -0500http://ask.sagemath.org/question/43369/Solve a differential equation using series expansions.http://ask.sagemath.org/question/42110/solve-a-differential-equation-using-series-expansions/Given an ODE such as $$y''+x^2y'+y=0$$ Is it possible to get sage to display the solution in the from (at least the first few terms of the expansion) $$y=a_o\left(c_0+c_1x+c_2x^2+\dots\right) + a_1\left(d_0+d_1x+d_2x^2+\dots\right) $$
[my attempts:](https://sagecell.sagemath.org/?z=eJxtzc2OwiAUhuG9iffwxS4EJcQ242rCFczW3UQNthCIHUgAK9z90C7mJ3FxNm-ec84kA9lkVlhiH9euTruh7-tVgYB-uD5Z78i2bCnJcx5U7YPVmhSWWUexR750u5_S0n0R4lBl84Ie_8G33ULRRKW-ImxCLx2iHyeFZGyEd4pzrFfR-CcZ1PzfzDfVYmphKAx2kqHWzNB7l4K9Xf2gxCk8loXGiM87DyYSypMsow-kypbhSKF9wB3WwZzRaGnHyH__mT_8MHP6DYNqVi0=&lang=sage)
EDIT:
I have made some progress, functional but it is not pretty.
[second attempt](https://sagecell.sagemath.org/?z=eJx1U9FumzAUfa_Uf7gqD7UTj2G6SlM6P-5pfpiUvUSIVU4ww8KxESYpaOq_z0BI6ZLKwrq6nHOur8_1UdToriUdaciP59h_lIiICHqHn25vOmCQH8yuUdag--4eo7ZPZ9LnM5XnqCMtiTEsIVq0v-PFOUfxsmMs8lhX2BeUyZ7GWTJIpz6Wlf8Lwfe2EiYDsbWHBn5aZRrwlDnQ8VBUlTQZclYfpZcis8o4idKwLhzCfYGSPfj9pVBaQvnt6-r2Bt7oA8vx5BNNhwP2BChZuaTTKR3HkwZnjo8n3ofusHXo76y1laAE2pWsCHSoZbLymYjADDFlKbxiyG0Ne1AGePqBZDwReBKnFwzNHqam9NjUNQ39pqEvNUAzPTS6ZglPyvRzLnaNrZXQqByxZY-thfkjkZYGcYzTJ_BmbFiy9gRvb_khzqMoQ-6wRxt8OpegK38lIlrRV-xzal9plXfIX_AmvoTSARr9Dx1c4VOwnoLNOaDnKPZRvwLvcJSyYT_pT05FvU2v-Crmvbsj6DwS0wquzP3ju5n_shimHgIn5d6BamAnDAyDC02hHFgjwxAm9fFZBKc2-7DoC8jzpENHQB1F7bMtgZ01Ta22zzaT7Fd9GMkFS8rxBYSN6LStkUf6-3yc2VqkEORCaRfOihczfNTj8T-h4CgD&lang=sage)userXFri, 20 Apr 2018 23:36:53 -0500http://ask.sagemath.org/question/42110/quick latex questionhttp://ask.sagemath.org/question/41644/quick-latex-question/ I am writing a worksheet for an ODE course and would like to display the following equation as $$y''+y=1$$ OR as $$\frac{d^2y}{dx^2}+y=1$$ However sage always reverts to partials. Is there an easy way to display these as ordinary differential equations?
thx
see sage cell [link text](https://sagecell.sagemath.org/?z=eJyrUFBQsFUoSyzSUK9Q17Tm5aoEC6SV5iWXZObnaahXqmtqVIAkXFwDgRIpmWlpGpU6FZoK2gqVCroKhra2BkDJ4vwckGQqkC5L1QAq1VGoBGkCymTkl4MENHm5Cooy80o0chJLUivAIpoAAaUhcw==&lang=sage)userXMon, 19 Mar 2018 13:25:22 -0500http://ask.sagemath.org/question/41644/Transparency option in streamline plothttp://ask.sagemath.org/question/37728/transparency-option-in-streamline-plot/In most graphics functions it's possible to set the transparency of the plot with the `alpha` option. In a streamline plot, e.g. take
sage: var('x, y')
sage: streamline_plot((y, -x**2 * y - x + y), (x, -3, 3), (y, -3, 3), density=1.4, color='black')
it doesn't work:
sage: var('x, y')
sage: streamline_plot((y, -x**2 * y - x + y), (x, -3, 3), (y, -3, 3), density=1.4, color='black', alpha=0.5)
...
WARNING: Ignoring option 'alpha'=0.500000000000000
...
Is there a workaround for this?
mforetsMon, 29 May 2017 08:20:26 -0500http://ask.sagemath.org/question/37728/solving a physic problem using sagehttp://ask.sagemath.org/question/10625/solving-a-physic-problem-using-sage/Hi, I'm new in this community.
I want to solve a physic problem which requires differential equation system solutions.
I don't know if my equations are correctly set. Any suggestion is good. My problem is described by this image: http://img805.imageshack.us/img805/7043/lllzm.png
I have two masses (1/3m the first, 2/3m the second) linked with a rope. The rope is free to slide around a nail (the big black point in the image). The image shows the starting condition: a man keeps the first mass stopped and so the rope is kept stretched by the second mass.
I search three functions describing the kinematics of two masses after the man will leave the fist mass: vertical movement of mass A y(t), vertical movement of mass B j(t), and horizontal movement of mass B x(t).
My Cartesian reference system is x-y system in the image.
I have to solve the following equations:
1. $-\frac{2}{3}mg+T=\frac{2}{3}m \frac{d^2y}{dt^2}$
2. $-\frac{1}{3}mg+S_y=\frac{1}{3}m\frac{d^2j}{dt^2}$
3. $S_x=\frac{1}{3}m\frac{d^2x}{dt^2}$
4. $|T|=\sqrt{S_x^2+S_y^2}$
5. $|y(t)|=\sqrt{x(t)^2+j(t)^2}$
From the forth and the fifth equations I obtain two equations, so I have 5 equations in 5 unknowns. They are:
1) T force sustaining the second mass
2) Sx x-component of force sustaining the first mass
3) Sy y-component of force sustaining the first mass
4) y(t) position of the second mass
5) x(t) x-position of the first mass
6) j(t) y-position of the first mass
I hope my explanation is clear.
How can I obtain my solutions using Sage?
Thank you very much!!pspFri, 18 Oct 2013 07:33:59 -0500http://ask.sagemath.org/question/10625/Symmetry-findinghttp://ask.sagemath.org/question/36844/symmetry-finding/On pp. 152-3 of Hydon's [*Symmetry Methods for Differential Equations*](https://www.kent.ac.uk/smsas/personal/ph282/sym.htm) (2000 ed.), he lists some computer packages for symmetry-finding. [This related Mathematica StackExchange question](https://mathematica.stackexchange.com/q/20438/47125) mentions the [SYM](https://www.math.upatras.gr/~spawn/SYM_html/guide/SYMOverview.html) Mathematica package and Maple's [DEtools/symgen](https://www.maplesoft.com/support/help/Maple/view.aspx?path=DEtools/symgen). Does SAGE have anything similar for doing symmetry-finding?GeremiaMon, 06 Mar 2017 12:27:37 -0600http://ask.sagemath.org/question/36844/Error: list' object is not callablehttp://ask.sagemath.org/question/36563/error-list-object-is-not-callable/I am trying to run the following code :
x = var('x')
y = function('y')(x)
de= x*x*diff(y,x,2)+x*diff(y,x)+((12)*x*x-4)*y == 0
val_test=desolve(de,y,contrib_ode=true)
a=val_test(4)
a.n()
But a.n() shows "Error: list' object is not callable".AeimiFri, 10 Feb 2017 23:33:43 -0600http://ask.sagemath.org/question/36563/Using desolve_odeint with external stimulushttp://ask.sagemath.org/question/36132/using-desolve_odeint-with-external-stimulus/I'm trying to numerically solve a two-variable system of differential equations with a periodic external driver (square wave). Right now, the code I have is:
from scipy.signal import square
def sqstim(t, amp=1):
return (amp/2)*float(square(t))+amp/2
var("V,w")
trange=srange(0,100,0.1)
def eqns(t):
return [1/100*(-w+(V^3-V)+sqstim(t)), V-0.2*w]
sol=desolve_odeint(eqns, times=trange, ics=[0.1,0.1], dvars=[V,w])
This gives `"TypeError: 'function' object is not iterable"`. What should I do instead?jaiaSat, 31 Dec 2016 15:54:09 -0600http://ask.sagemath.org/question/36132/Pattern matching in differential equationshttp://ask.sagemath.org/question/10452/pattern-matching-in-differential-equations/We have the following functional equation:
f(x*y)+f(x*(1-y))+f((1-x)*y)+f((1-x)*(1-y)) == f(x)+f(1-x)+f(y)+f(1-y)
Differentiating by x then y we get a second order differential equation which contains the pattern g(t)=f'(t)+tf''(t) three times with different expressions for t.
How can I change the occurences of the patterns to the appropriate g(.)?
czsanSun, 18 Aug 2013 04:04:58 -0500http://ask.sagemath.org/question/10452/Syntax for numerically solve differential equationhttp://ask.sagemath.org/question/34096/syntax-for-numerically-solve-differential-equation/ I am having trouble using the desolve_rk4 function and I assume I am messing up the syntax somewhere. I tried
sage: t = var('t')
sage: g = function('g')(t)
sage: beta(x) = (8*x^3*(1 + x))/(-1/8 + x)
sage: eq = diff(g, t) - 2*pi*I*beta(x=g)/diff(beta, x)(x=g)
sage: desolve_rk4(eq, g, ics=[0, 0.139202795060905 + 0.0642260199336103*I], end_points=[1])
but it just gives me back a constant solution:
[[0, -0.363218465306351],
[0.1, -0.363218465306351],
[0.2, -0.363218465306351],
[0.3, -0.363218465306351],
[0.4, -0.363218465306351],
[0.5, -0.363218465306351],
[0.6000000000000001, -0.363218465306351],
[0.7000000000000001, -0.363218465306351],
[0.8, -0.363218465306351],
[0.9, -0.363218465306351],
[1.0, -0.363218465306351]]
The only other thing I have seen is that Sage maybe can't handle complex functions.
If it helps, the corresponding Mathematica code is
beta[x_] := 8 x^3 (1 + x)/(-1/8 + x)
NDSolve[{g'[t] == 2 Pi I beta[g[t]]/beta'[g[t]], g[0] == -0.363218465306351}, g, {t, 0, 1}]
jaebondFri, 15 Jul 2016 13:16:23 -0500http://ask.sagemath.org/question/34096/Solving system of ODEshttp://ask.sagemath.org/question/27248/solving-system-of-odes/I'm trying to solve a system of 4 ODEs using Sage. I'm able to get a solution, but it's not at all what I expect. Here are the 4 equations:
$$\frac{dw}{dt} = \frac{Q_x}{V} x(t) - \frac{Q_y}{V}\eta w(t) - \frac{Q_z}{V} \eta w(t)$$
$$\frac{dx}{dt} = \frac{Q_y}{V} \gamma y(t) + \frac{Q_z}{V} \zeta z(t) - \frac{Q_x}{V}x(t)$$
$$\frac{dy}{dt} = \frac{Q_y}{V} \eta w(t) - \frac{Q_y}{V} \gamma y(t) $$
$$\frac{dz}{dt} = \frac{Q_z}{V} \eta w(t) - \frac{Q_z}{V} \zeta z(t)$$
Here is how I've coded the problem in Sage.
eta, zeta, gamma = var('eta, zeta, gamma')
m_0, Qx, Qy, Qz, V = var('m_0, Q_x, Q_y, Q_z, V')
t = var('time')
w=function('w', t)
x=function('x', t)
y=function('y', t)
z=function('z', t)
dw = diff(w,t) - Qx/V*x + Qy/V*eta*w + Qz/V*eta*w == 0
dx = diff(x,t) - Qy/V*gamma*y - Qz/V*zeta*z + Qx/V*x == 0
dy = diff(y,t) - Qy/V*eta*w + Qy/V*gamma*y == 0
dz = diff(z,t) - Qz/V*eta*w + Qz/V*zeta*z == 0
eqs = [dw, dx, dy, dz]
vars = [w, x, y, z]
ics = [0, 0, m_0, 0, 0]
ans = desolve_system(eqs, vars, ics, ivar=t)
print ans
Where m_0 is an initial mass that exists in the x volume at time 0. All other volumes are empty at time 0.
I'm expecting relatively simple exponential equations as the solution, but what I get in return is some significantly more complex.
$$
w\left(\mathit{time}\right) = \frac{V \gamma m_{0} \zeta}{{\left(\eta \gamma + {\left(\eta \gamma + \eta\right)} \zeta\right)} V + \gamma \zeta} + \mathcal{L}^{-1}\left(-\frac{V^{4} g_{2562}^{2} \gamma m_{0} \zeta + {\left(Q_{x} Q_{z} - Q_{z}^{2}\right)} V^{2} \gamma^{2} m_{0} \zeta^{2} - {\left(Q_{x} Q_{z} \eta m_{0} \zeta^{2} + {\left(Q_{x}^{2} - Q_{x} Q_{z}\right)} \eta \gamma^{2} m_{0}\right)} V^{3} - {\left({\left(Q_{x} \eta \gamma m_{0} + Q_{x} \eta m_{0} \zeta\right)} V^{4} - {\left({\left(Q_{x} - Q_{z}\right)} \gamma^{2} m_{0} \zeta + Q_{z} \gamma m_{0} \zeta^{2}\right)} V^{3}\right)} g_{2562}}{{\left(V^{3} g_{2562}^{3} + {\left(Q_{x} V^{3} \eta + {\left({\left(Q_{x} - Q_{z}\right)} \gamma + Q_{z} \zeta + Q_{x}\right)} V^{2}\right)} g_{2562}^{2} + {\left(Q_{x}^{2} Q_{z} - Q_{x} Q_{z}^{2}\right)} \gamma \zeta + {\left({\left(Q_{x}^{2} Q_{z} - Q_{x} Q_{z}^{2}\right)} \eta \gamma + {\left({\left(Q_{x}^{2} Q_{z} - Q_{x} Q_{z}^{2}\right)} \eta \gamma + {\left(Q_{x}^{2} Q_{z} - Q_{x} Q_{z}^{2}\right)} \eta\right)} \zeta\right)} V + {\left({\left(Q_{x} Q_{z} \eta \zeta + Q_{x}^{2} \eta + {\left(Q_{x}^{2} - Q_{x} Q_{z}\right)} \eta \gamma\right)} V^{2} + {\left({\left(Q_{x}^{2} - Q_{x} Q_{z}\right)} \gamma + {\left(Q_{x} Q_{z} + {\left(Q_{x} Q_{z} - Q_{z}^{2}\right)} \gamma\right)} \zeta\right)} V\right)} g_{2562}\right)} {\left({\left(\eta \gamma + {\left(\eta \gamma + \eta\right)} \zeta\right)} V + \gamma \zeta\right)}}, g_{2562}, \mathit{time}\right)
\end{equation}
\begin{equation}
x\left(\mathit{time}\right) = \frac{V \eta \gamma m_{0} \zeta}{{\left(\eta \gamma + {\left(\eta \gamma + \eta\right)} \zeta\right)} V + \gamma \zeta} + \mathcal{L}^{-1}\left(\frac{{\left(Q_{x} Q_{z} - Q_{z}^{2}\right)} V \gamma^{2} m_{0} \zeta^{2} + {\left(Q_{x} Q_{z} \eta^{2} m_{0} \zeta^{2} + {\left(Q_{x}^{2} - Q_{x} Q_{z}\right)} \eta^{2} \gamma^{2} m_{0}\right)} V^{3} + {\left({\left(Q_{x} Q_{z} - Q_{z}^{2}\right)} \eta \gamma^{2} m_{0} \zeta + {\left(Q_{x} Q_{z} - Q_{z}^{2}\right)} \eta \gamma m_{0} \zeta^{2}\right)} V^{2} + {\left(V^{3} \gamma m_{0} \zeta + {\left(\eta \gamma m_{0} + \eta m_{0} \zeta\right)} V^{4}\right)} g_{2562}^{2} + {\left({\left(Q_{x} \eta^{2} \gamma m_{0} + Q_{x} \eta^{2} m_{0} \zeta\right)} V^{4} + {\left({\left(Q_{x} - Q_{z}\right)} \eta \gamma^{2} m_{0} + Q_{x} \eta \gamma m_{0} \zeta + Q_{z} \eta m_{0} \zeta^{2}\right)} V^{3} + {\left({\left(Q_{x} - Q_{z}\right)} \gamma^{2} m_{0} \zeta + Q_{z} \gamma m_{0} \zeta^{2}\right)} V^{2}\right)} g_{2562}}{{\left(V^{3} g_{2562}^{3} + {\left(Q_{x} V^{3} \eta + {\left({\left(Q_{x} - Q_{z}\right)} \gamma + Q_{z} \zeta + Q_{x}\right)} V^{2}\right)} g_{2562}^{2} + {\left(Q_{x}^{2} Q_{z} - Q_{x} Q_{z}^{2}\right)} \gamma \zeta + {\left({\left(Q_{x}^{2} Q_{z} - Q_{x} Q_{z}^{2}\right)} \eta \gamma + {\left({\left(Q_{x}^{2} Q_{z} - Q_{x} Q_{z}^{2}\right)} \eta \gamma + {\left(Q_{x}^{2} Q_{z} - Q_{x} Q_{z}^{2}\right)} \eta\right)} \zeta\right)} V + {\left({\left(Q_{x} Q_{z} \eta \zeta + Q_{x}^{2} \eta + {\left(Q_{x}^{2} - Q_{x} Q_{z}\right)} \eta \gamma\right)} V^{2} + {\left({\left(Q_{x}^{2} - Q_{x} Q_{z}\right)} \gamma + {\left(Q_{x} Q_{z} + {\left(Q_{x} Q_{z} - Q_{z}^{2}\right)} \gamma\right)} \zeta\right)} V\right)} g_{2562}\right)} {\left({\left(\eta \gamma + {\left(\eta \gamma + \eta\right)} \zeta\right)} V + \gamma \zeta\right)}}, g_{2562}, \mathit{time}\right)
\end{equation}
\begin{equation}
y\left(\mathit{time}\right) = \frac{V \eta m_{0} \zeta}{{\left(\eta \gamma + {\left(\eta \gamma + \eta\right)} \zeta\right)} V + \gamma \zeta} + \mathcal{L}^{-1}\left(-\frac{V^{4} \eta g_{2562}^{2} m_{0} \zeta + {\left(Q_{x} Q_{z} \eta + {\left(Q_{x} Q_{z} - Q_{z}^{2}\right)} \eta \gamma\right)} V^{2} m_{0} \zeta^{2} + {\left(Q_{x} Q_{z} \eta^{2} m_{0} \zeta^{2} + Q_{x} Q_{z} \eta^{2} m_{0} \zeta - {\left(Q_{x}^{2} - Q_{x} Q_{z}\right)} \eta^{2} \gamma m_{0}\right)} V^{3} + {\left(Q_{x} V^{4} \eta^{2} m_{0} \zeta + {\left(Q_{z} \eta m_{0} \zeta^{2} + {\left({\left(Q_{x} - Q_{z}\right)} \eta \gamma + Q_{x} \eta\right)} m_{0} \zeta\right)} V^{3}\right)} g_{2562}}{{\left(V^{3} g_{2562}^{3} + {\left(Q_{x} V^{3} \eta + {\left({\left(Q_{x} - Q_{z}\right)} \gamma + Q_{z} \zeta + Q_{x}\right)} V^{2}\right)} g_{2562}^{2} + {\left(Q_{x}^{2} Q_{z} - Q_{x} Q_{z}^{2}\right)} \gamma \zeta + {\left({\left(Q_{x}^{2} Q_{z} - Q_{x} Q_{z}^{2}\right)} \eta \gamma + {\left({\left(Q_{x}^{2} Q_{z} - Q_{x} Q_{z}^{2}\right)} \eta \gamma + {\left(Q_{x}^{2} Q_{z} - Q_{x} Q_{z}^{2}\right)} \eta\right)} \zeta\right)} V + {\left({\left(Q_{x} Q_{z} \eta \zeta + Q_{x}^{2} \eta + {\left(Q_{x}^{2} - Q_{x} Q_{z}\right)} \eta \gamma\right)} V^{2} + {\left({\left(Q_{x}^{2} - Q_{x} Q_{z}\right)} \gamma + {\left(Q_{x} Q_{z} + {\left(Q_{x} Q_{z} - Q_{z}^{2}\right)} \gamma\right)} \zeta\right)} V\right)} g_{2562}\right)} {\left({\left(\eta \gamma + {\left(\eta \gamma + \eta\right)} \zeta\right)} V + \gamma \zeta\right)}}, g_{2562}, \mathit{time}\right)
\end{equation}
\begin{equation}
z\left(\mathit{time}\right) = \frac{V \eta \gamma m_{0}}{{\left(\eta \gamma + {\left(\eta \gamma + \eta\right)} \zeta\right)} V + \gamma \zeta} + \mathcal{L}^{-1}\left(-\frac{V^{4} \eta g_{2562}^{2} \gamma m_{0} - {\left(Q_{x} Q_{z} \eta^{2} m_{0} \zeta - {\left({\left(Q_{x}^{2} - Q_{x} Q_{z}\right)} \eta^{2} \gamma^{2} + {\left(Q_{x}^{2} - Q_{x} Q_{z}\right)} \eta^{2} \gamma\right)} m_{0}\right)} V^{3} + {\left({\left(Q_{x} Q_{z} - Q_{z}^{2}\right)} \eta \gamma^{2} m_{0} \zeta + {\left(Q_{x}^{2} - Q_{x} Q_{z}\right)} \eta \gamma^{2} m_{0}\right)} V^{2} + {\left(Q_{x} V^{4} \eta^{2} \gamma m_{0} + {\left(Q_{z} \eta \gamma m_{0} \zeta + {\left({\left(Q_{x} - Q_{z}\right)} \eta \gamma^{2} + Q_{x} \eta \gamma\right)} m_{0}\right)} V^{3}\right)} g_{2562}}{{\left(V^{3} g_{2562}^{3} + {\left(Q_{x} V^{3} \eta + {\left({\left(Q_{x} - Q_{z}\right)} \gamma + Q_{z} \zeta + Q_{x}\right)} V^{2}\right)} g_{2562}^{2} + {\left(Q_{x}^{2} Q_{z} - Q_{x} Q_{z}^{2}\right)} \gamma \zeta + {\left({\left(Q_{x}^{2} Q_{z} - Q_{x} Q_{z}^{2}\right)} \eta \gamma + {\left({\left(Q_{x}^{2} Q_{z} - Q_{x} Q_{z}^{2}\right)} \eta \gamma + {\left(Q_{x}^{2} Q_{z} - Q_{x} Q_{z}^{2}\right)} \eta\right)} \zeta\right)} V + {\left({\left(Q_{x} Q_{z} \eta \zeta + Q_{x}^{2} \eta + {\left(Q_{x}^{2} - Q_{x} Q_{z}\right)} \eta \gamma\right)} V^{2} + {\left({\left(Q_{x}^{2} - Q_{x} Q_{z}\right)} \gamma + {\left(Q_{x} Q_{z} + {\left(Q_{x} Q_{z} - Q_{z}^{2}\right)} \gamma\right)} \zeta\right)} V\right)} g_{2562}\right)} {\left({\left(\eta \gamma + {\left(\eta \gamma + \eta\right)} \zeta\right)} V + \gamma \zeta\right)}}, g_{2562}, \mathit{time}\right)
$$
Can anyone explain what I've done wrong?spontarelliamThu, 02 Jul 2015 13:52:10 -0500http://ask.sagemath.org/question/27248/Problems in solving eigenvalue equations with differential operatorshttp://ask.sagemath.org/question/29398/problems-in-solving-eigenvalue-equations-with-differential-operators/Let's say we have a linear differential operator $\hat A = \sum_k a_k \frac {d^k} {dx^k}$ and we want to solve the eigenvalue equation $\hat A f(x) = a f(x)$ which is an ODE that we put in sage and outputs a solution $g(x)$. We now want to substitute $f(x)$ with $g(x)$ and simplify our expression in order to extract the eigenvalues of $\hat A$. The problem here is that I am not able to substitute the derivatives with my solution neither via `eqn.substitute_expression(f(x) == g(x))` nor via `eqn.substitute_function(f(x),g(x))` because `D[0]f(x)`, ... `D[0,0,...,0]f(x)` remain unchanged.
papachristoumariosThu, 03 Sep 2015 16:50:54 -0500http://ask.sagemath.org/question/29398/SageManifold solving differential equations (Schwarzschild metric)http://ask.sagemath.org/question/28774/sagemanifold-solving-differential-equations-schwarzschild-metric/ Hi
I've been following Schwarzschild derivation here https://en.wikipedia.org/wiki/Deriving_the_Schwarzschild_solution#Using_the_field_equations_to_find_A.28r.29_and_B.28r.29 and I've come to the point where I have the three differential equations listed in the above link.
I tried to solve it with
> desolve_system
as follows:
https://www.dropbox.com/s/1uihjhzm2k1wl4a/2015-08-07-133804.pdf?dl=0
But it does not seem to be working. Am I missing something? It's a system of nonlinear differential equationsOttoSat, 08 Aug 2015 01:08:26 -0500http://ask.sagemath.org/question/28774/Solving this DE containing an integralhttp://ask.sagemath.org/question/28643/solving-this-de-containing-an-integral/I am trying to solve
$$ 0 = - \partial_a F(a)-e(a)F(a) + \int_0^{a} e(a')\partial_a F(a') d a' + n $$
optimally, without specifying e(a). But if necessary (as I guess), e(a) = k1*exp(k2*a). Here's my code:
var('a b k_1 k_2 n')
e(b) = k_1*exp(b*k_2)
F = function('F', a)
g(b) = e(b)*diff(F,b,1)
de = -diff(F, a, 1) - e(a)*F(a) + g(b).integral(b, 0, a) + n
y = desolve(de, F, ivar=a); y
And the output is
(_C + n*Ei(k_1*e^(a*k_2)/k_2)/k_2)*e^(-k_1*e^(a*k_2)/k_2)
However, I believe that something is wrong with my integral, it is not doing what I expect it to do. For example, if I change the integration boundary to 1:
de = -diff(F, a, 1) - e(a)*F(a) + g(b).integral(b, 0, a) + n
I will still get the same result:
(_C + n*Ei(k_1*e^(a*k_2)/k_2)/k_2)*e^(-k_1*e^(a*k_2)/k_2)
Where am I going wrong?FooBarSun, 19 Jul 2015 07:13:44 -0500http://ask.sagemath.org/question/28643/Solve DE with secondary functionhttp://ask.sagemath.org/question/26820/solve-de-with-secondary-function/ I want to solve
$-F''(a) + g(a)( F'(a) - F(a)$
where $g(a)$ is logistic.
My code is
a = var('a')
g(a)= 1 /(1 + e^(-(a)))
F = function('F', a)
de = -diff(F,a, 2) + g(a)*(diff(F,a,1) - F(a))
y = desolve(de, F); y
But I get an error (at the end of the post). If I remove `e(a)` from the `de`-term, I get a solution. What do I have to take into account when I add the second function?
TypeError: no canonical coercion from <class 'sage.symbolic.function_factory.NewSymbolicFunction'> to Symbolic RingFooBarWed, 13 May 2015 09:41:29 -0500http://ask.sagemath.org/question/26820/Compact output of solution of DEhttp://ask.sagemath.org/question/26316/compact-output-of-solution-of-de/ When I'm trying to solve DE:
t = var('t')
y = function('y', t)
de = t*(y^2)*diff(y,t) + y^3 == 1
sol = desolve(de,[y,t], [1,2])
the output is pretty ugly:
-1/3*log(y(t)^2 + y(t) + 1) - 1/3*log(y(t) - 1) == -1/3*log(7) + log(t)
When I'm solving this in matlab:
clear;
syms y(t)
y(t) = dsolve(t*(y^2)*diff(y,t) + y^3 == 1, y(1) ==2)
The output looks much better:
y(t) = (exp(log(7) - 3*log(t)) + 1)^(1/3)
Can I see output in sage looking similiar to this from matlab? Simplify(sol) dosen't work. Maybe I've made mistake somewhere, but I can't determine without knowing the form y(t) from sage.
And btw, typing:
t*(y^2)*y'+ y^3 = 1, y(1) = 2
into wolframalpha.com results yet another solution. I'm lost...PhotonTue, 24 Mar 2015 15:51:52 -0500http://ask.sagemath.org/question/26316/Phase portraits of 2-dimensional systemshttp://ask.sagemath.org/question/9423/phase-portraits-of-2-dimensional-systems/I'm trying to plot solutions to two dimensional ordinary differential equations. I found that Sage makes it easy to plot a vector field and, using `ode_solver()`, I can plot solutions on top of the vector field by specifying an initial condition `y_0` and some range of time to run (`t_span`).
However, this method I'm using seems to be quite ad hoc, as I have to choose the right initial conditions and time span / know a lot about my system in order to plot a nice picture. Let's make this more concrete:
Say I want to draw a nice phase portrait for
$\dot{x} = -y$
$\dot{y} = -x$
First I generate the vector field:
var('x y')
VF=plot_vector_field([-y,-x],[x,-2,2],[y,-2,2])
Then I use `ode_solver()` to plot solutions with initial conditions going around the unit circle:
T = ode_solver()
T.function=lambda t,y:[-y[1],-y[0]]
solutions = []
c = circle((0,0), 1, rgbcolor=(1,0,1))
for k in range(0,8):
T.ode_solve(y_0=[cos(k*pi/4),sin(k*pi,t_span=[0,1],num_points=100)
solutions.append(line([p[1] for p in T.solution]))
This generates the following picture:
![](http://oi46.tinypic.com/27y4v80.jpg)
But if I change run the system for one more unit of time (set `t_span=[0,2]`), the picture gets skewed:
![](http://oi49.tinypic.com/2mfj1i8.jpg)
This makes sense, of course, because the magnitude of the vectors along $y=-x$ get big as you get further away from the origin. Similarly, the trajectory along $y=x$ has trouble getting to the origin because the magnitude of those vectors get very small. Which all makes me think there's a better way to do this. Thoughts?dxvxdSun, 14 Oct 2012 08:02:40 -0500http://ask.sagemath.org/question/9423/Using lists to create function to automate Runge-Kutta Methodhttp://ask.sagemath.org/question/25698/using-lists-to-create-function-to-automate-runge-kutta-method/I'm in a differential equations class in which we're expected to use computers to solve some of the problems. My professor uses Maple, but I'm choosing to use Sage because it won't cost me hundreds of dollars. I have a very limited CS background.
One part of our homework this week is to estimate the value of y(x_n) given dy/dx = f(x,y) and y(x_0) = y_0 using the Runge-Kutta Method. I know this shouldn't be that hard to do in Sage, but I'm having a hard time figuring this out. Any pointers? r.jay.hutchinsonTue, 03 Feb 2015 16:33:20 -0600http://ask.sagemath.org/question/25698/Design of classes for dynamical systemshttp://ask.sagemath.org/question/25403/design-of-classes-for-dynamical-systems/ I have an ongoing project using Sage to process some dynamical systems (ODEs, that is), and along the way I've been developing some classes to generalize the process of storing an ODE system and working with it. I'd be interested in contributing the work to Sage if it's helpful, and in understanding what's a "sage-like" way to do it.
So here I'm wondering about what would be worth trying to contribute to Sage, and asking whether people would do the design differently.
I have a central class called ODEsystem, which is built around a dictionary encoding the flow function, e.g. if the system is
$$dx/dt = -ax-by$$
$$dy/dt = -cx-dy$$
the flow is a dictionary `{ x:-ax-by, y:-cx-dy }` whose keys are variables and values are symbolic expressions. The object has methods to write itself out in LaTeX form, to integrate the dynamics from initial conditions, to find its equilibria and their stability, to make bifurcation diagrams under certain conditions, and a few other things. My intention is to allow room for a class hierarchy of models to develop, including difference equations, stochastic dynamics, PDEs, etc. There's also room to create subclasses and functions that derive specific models from general ones, simple models as limit cases of complex models, etc.
Another class that may be of use to others is Bindings, which binds variables to specific values in various expressions. For example, if I write
var('a, b')
generic_system = ODESystem( { x:a-bx }, vars=[x] )
specific_bindings = Bindings( a=1, b=2 )
latex( generic_system.bind( specific_bindings ) )
I'll get $dx/dt = 1-2x$. This is a wrapper for `substitute()` and `substitute_function()`, of course, but adds convenience and flexibility. (Note ODEsystem uses 'vars' to define the order of state variables. It's not logically necessary but it's useful in multiple dimensions.)
Code for those two classes is at https://github.com/worden-lee/SageDynamics/blob/master/SageDynamics/dynamicalsystems.py.step, and there's some more stuff in other files there. Example output is at http://lalashan.mcmaster.ca/theobio/worden/index.php/SageDynamics/Demo.
Questions:
* Could something like this become a useful part of Sage
* Would people want it to be designed differently
* Different situations call for different choices of integration routines. One strategy is to provide different ODEsystem classes that use different numerical integrators, and another one is to uncouple the integration from the class that represents the system. I'm not sure what would be better.
* Would it be good to represent the ODE as an element of a ring of vector fields on an appropriate manifold, with some extra methods added? It seems like that would provide some power, but I don't know the system well enough to go there right now.
* Should I send this to sage-devel or something instead of posting herewonderMon, 05 Jan 2015 14:45:15 -0600http://ask.sagemath.org/question/25403/solving a system of DEs numerically and plotting the solutionhttp://ask.sagemath.org/question/24766/solving-a-system-of-des-numerically-and-plotting-the-solution/Hi everyone. I want to hand Sage a system of nonlinear DEs with initial values, and I'd like a plot of the solutions. Is there a best way to do this? For example, to solve this system: ds/dt=-si, di/st=si-2i, dr/dt=2i, s(0)=1, i(0)= 0.00000127, r(0)=0, I did this in Sage:
sage: (s,i,r)=var('s,i,r')
sage: des=[-1*s*i, s*i-2*i, i]
sage: desolve_system_rk4(des, [s,i,r], ics=[0, 1, 0.00000127, 0], ivar=t, end_points=20)
Sage gave me a list of points, but I couldn't figure out a way to plot s, i, and r using the points it gave me. Thanks!
strangelove1661Mon, 03 Nov 2014 14:32:05 -0600http://ask.sagemath.org/question/24766/Differential Equations and Error of Estimationshttp://ask.sagemath.org/question/24330/differential-equations-and-error-of-estimations/ So I'm tasked with the following two questions:
1. For the initial value problem yβ²=y sin(x), y(0)=1, first find the exact solution. Then make log-log plots of the error versus n (the number of steps) at x=pi for Euler's method and the fourth-order Runge-Kutta method.
2. Repeat the above exercise for the IVP yβ²=y+sin(x), y(0)=β0.48.
So we have a defined function for funding the Euler estimation:
def EulerMethod(xstart, ystart, xfinish, nsteps, f):
'''
Returns a list of x and y values for the initial value problem
y' = f, y(xstart) = ystart, up to x=xfinish, using n steps of
Euler's Method.
EXAMPLE:
var('x,y')
f(x,y) = -y + cos(x)
EulerMethod(0,1,2,4,f)
[(0, 1), (0.500000000000000, 1), (1.00000000000000, 0.938791280945),
(1.50000000000000, 0.739546793407), (2.00000000000000, 0.405141997537)]
#If you just want to get the last value, you could do:
sol = EulerMethod(0,1,2,4,f)
sol[-1]
(2.00000000000000, 0.405141997537)
'''
sol = [ystart]
xvals = [xstart]
h = N((xfinish-xstart)/nsteps)
for step in range(nsteps):
sol.append(sol[-1] + h*f(xvals[-1],sol[-1]))
xvals.append(xvals[-1] + h)
return zip(xvals,sol)
As written by the teacher. We are given that and no function for the RK4 so I assume we are to use the built in Sage RK4 estimation.
Anyways I have no idea what to do. I'm absolutely lost. Can anyone please help me? Massive thanks in advance.HankSpankMon, 29 Sep 2014 17:45:08 -0500http://ask.sagemath.org/question/24330/How to graph this function in Sage plot?http://ask.sagemath.org/question/24171/how-to-graph-this-function-in-sage-plot/ This function:
clogxβxd=βalogy+by+y0
Thank you.jukhamilWed, 17 Sep 2014 17:02:14 -0500http://ask.sagemath.org/question/24171/