ASKSAGE: Sage Q&A Forum - RSS feedhttps://ask.sagemath.org/questions/Q&A Forum for SageenCopyright Sage, 2010. Some rights reserved under creative commons license.Mon, 19 Feb 2024 16:46:11 +0100Comparison between SageMath ODE solvers and Scipy ODE solvershttps://ask.sagemath.org/question/76091/comparison-between-sagemath-ode-solvers-and-scipy-ode-solvers/ I'm pretty new to SageMath and I'm currently trying to work on solving geodesic equations in well-defined spacetime manifolds. I found the documentation of SageMath pretty extensive and broad and it gave me a lot of confidence in the software. But I'm still hesitating to use Sage ODE solvers because I know the capabilities of `RK45` from `scipy.integrate` as I have use that many times. I just want to know on a broad scale, are SageMath solvers faster than Scipy solvers and how reliable and stable are the solutions in general ?
Thanks in advance for the replies! AKGSageMon, 19 Feb 2024 16:46:11 +0100https://ask.sagemath.org/question/76091/Error With Odeinthttps://ask.sagemath.org/question/69589/error-with-odeint/ v=3
t,x,y,dx,dy=var('t,x,y,dx,dy')
r=sqrt((x-v*t)^2+y^2)
f=(tanh(r+4)-tanh(r-4))/(2*tanh(4))
ft=derivative(f,t)
fx=derivative(f,x)
fy=derivative(f,y)
t_tt=f^2*fx*v^3
t_tx=-f*fx*v^2
t_xx=v*fx
t_ty=-(f*fy*v^2)/2
t_xy=fy*v/2
x_tt=f^3*fx*v^4-f*fx*v^2-v*ft
x_tx=-t_tt
x_xx=-t_tx
x_ty=-(f^2*fy*v^3+v*fy)/2
x_xy=(f*fy*v^2)/2
y_tt=-f*fy*v^2
y_tx=v*fy/2
g=-x_tt-x_tx*dx-x_xx*dx^2-x_ty*dy-x_xy*dx*dy+dx*(t_tt+t_tx*dx+t_xx*dx^2+t_ty*dy+t_xy*dx*dy)
h=-y_tt-y_tx*dx+dy*(t_tt+t_tx*dx+t_xx*dx^2+t_ty*dy+t_xy*dx*dy)
z=[dx,dy,g,h]
sol=desolve_odeint(z,[0.1,0.1,0,1],srange(0,10,0.1),[x,y,dx,dy],ivar=t)
p=line(zip(sol[:,0],sol[:,1]))
p.show()
This code gives me an error saying "excess work done," so I'm assuming that the derivatives blew up somewhere, but I don't know why it says this for these initial conditions and not for others. Now, odeint has parameters hmin, hmax, rtol, and atol, but I don't know if there is a way to fix a problem using these.Jack ZuffanteThu, 29 Jun 2023 00:16:51 +0200https://ask.sagemath.org/question/69589/Running desolve_odeint with full_output=1https://ask.sagemath.org/question/69413/running-desolve_odeint-with-full_output1/Sometimes running desolve_odeint gives a message that excess work was done and it recommends running with full_output=1. However, I haven't gotten that to work. Is it something else, like verbose=True?Jack ZuffanteFri, 23 Jun 2023 23:44:03 +0200https://ask.sagemath.org/question/69413/Solutions to desolve_odeint seem largehttps://ask.sagemath.org/question/69391/solutions-to-desolve_odeint-seem-large/I'm trying to figure out what, if anything, is going wrong with this code here. It says in the SageMath documentation that for desolve_odeint, the order of dvars must be the same as that of des. In my case, that means dp/dt=dp, dth/dt=dth,
d(dp)/t=g, and d(dth)/t=h. But it's spitting out values that are much larger than what I would expect for t-values between 0 and 10.
t,p,th,dp,dth=var('t,p,th,dp,dth')
g=5*p*dth^2+4*t*dth^2*dp
h=-4*t/(5*p^2+4*t^2)*dth-5*p/(5*p^2+4*t^2)*dp*dth+4*t*dth^3
f=[dp,dth,g,h]
sol=desolve_odeint(f,[5,0,1,1],srange(0,10),[p,th,dp,dth],ivar=t)
sol[2]Jack ZuffanteFri, 23 Jun 2023 04:35:56 +0200https://ask.sagemath.org/question/69391/System of Second-Order Differential Equations Part 2https://ask.sagemath.org/question/69150/system-of-second-order-differential-equations-part-2/ This is one way I know of to solve a system of differential equations.
T = ode_solver()
f = lambda t, y: [y[2], y[3], -y[1] - 1, -y[0]]
T.function = f
T.ode_solve(y_0=[1, 1, 0, 0], t_span=[0, 20], num_points=1000)
Just one thing: if the -y[1]-1 and y[0] are replaced with much larger functions of y[0], y[1]. y[2], and y[3], how can I declare those functions before defining T?Jack ZuffanteWed, 14 Jun 2023 03:31:57 +0200https://ask.sagemath.org/question/69150/desolve_system_rk4 vs desolve_odeinthttps://ask.sagemath.org/question/36199/desolve_system_rk4-vs-desolve_odeint/I am getting different results (a result vs no result) from desolve_system_rk4 and desolve_odeint for the same system. What am I missing? My code is the following:
reset()
from sage.calculus.desolvers import desolve_system_rk4,desolve_odeint
var('eta,r,ph');
drdt=-sqrt(2.0*r + 200.0 / r - 100.0)*(r - 2.0) / (r^2)
dphdt=10.0*(r - 2.0) / (r^3)
sol1=desolve_system_rk4([drdt,dphdt],[r,ph],ics=[0,4.0,0.3],ivar=eta,end_points=1000,step=0.01)
sol2=desolve_odeint([drdt,dphdt],[4.0,0.3],srange(0,1000,0.01),[r,ph])
D=[[j*cos(k),j*sin(k)] for i,j,k in sol1]
p1=list_plot(D,aspect_ratio=1)
p2=line(zip(sol2[:,0]*cos(sol2[:,1]),sol2[:,0]*sin(sol2[:,1])))
show(graphics_array([p1,p2]))implicitnoneFri, 06 Jan 2017 14:29:26 +0100https://ask.sagemath.org/question/36199/Using desolve_odeint with external stimulushttps://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 22:54:09 +0100https://ask.sagemath.org/question/36132/