Numerical Integration of System + x/y Plot

asked 2023-06-04 06:45:09 +0100

Jack Zuffante gravatar image

updated 2023-06-04 06:45:31 +0100

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 flag offensive close merge delete

Comments

I'll post my specific code tomorrow.

Jack Zuffante gravatar imageJack Zuffante ( 2023-06-04 08:22:11 +0100 )edit

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)

Jack Zuffante gravatar imageJack Zuffante ( 2023-06-04 18:49:23 +0100 )edit

Yeah, I'm not sure how to make it work

Jack Zuffante gravatar imageJack Zuffante ( 2023-06-04 20:02:21 +0100 )edit

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[2],y[3],-y[1]-1,-y[0]]

This is because the indices start with zero in SageMath. So y1 is y[0], y2 is y[1], 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[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)
f = T.interpolate_solution()
plot(f,0,5).show()

I hope that helps.

tolga gravatar imagetolga ( 2023-06-05 08:19:36 +0100 )edit