# Numerical Integration of System + x/y Plot

How do I numerically integrate such a system for a certain range of t and then make an x/y plot of the curve?

x''(t)=A(t,x,y,x',y')

y''(t)=B(t,x,y,x',y')

A and B are functions t, x, y, dx/dt, and dy/dt.

edit retag close merge delete

Well first off all theres the fact that just typing the double quotes there doesn't work. But my main problem is that in the example achrzesz posted, the dependent variable (which it looks like it is also t in his case) does not appear in the equations, so I guess I need some way to define t as a variable and x and y as functions of t in a way that would work with the example. My instinct is to do (using simpler equations for now):

t=var('t')

x=function('x')(t)

y=function('y')(t)

u=diff(x,t)

v=diff(y,t)

diff(u,t)=y*u

diff(v,t)=x*v

T=ode_solver()

f = lambda t,g:[u,v,diff(u,t),diff(v,t)]

T.function=f

T.ode_solve(g_0=[t0,x0,y0,dx0,dy0],t_span=[0,20],num_points=1000)

I would just ...(more)

The first part of that comment is showing the procedure that you need to do on the paper (reducing the second-order system to a system of the first-order).

The following part might be frustrating. After having your first-order equations as

y1'=y3 y2'=y4 y3'=-y2-1 y4'=-y1

you write in the code:

f = lambda t,y:[y,y,-y-1,-y]


This is because the indices start with zero in SageMath. So y1 is y, y2 is y, etc. (You think of y1, y2, etc as element of the array y)

The rest is the code. You will run this:

T=ode_solver()
f = lambda t,y:[y,y,-y-1,-y]
T.function=f
T.ode_solve(y_0=[1,1,0,0],t_span=[0,20],num_points=1000)
f = T.interpolate_solution()
plot(f,0,5).show()


I hope that helps.