numerical estimate of 2nd Order Diff Eq's in CAS

[Wes Loewer]

Is there a built-in method in CAS to numerically solve a 2nd Order Diff Eq, similar to the way odesolve() and plotode() do for a 1st Order Diff Eq?

[Brian Zilli]

I don't know if this qualifies as built-in since it does require some pre-processing, but this is one way to do it:

Suppose we have a second-order equation x''=f(t,x,x'). If we let x0=x and x1=x', we may transform this into a system of first-order equations:

x0' = x1
x1' = f(t,x0,x1)

In Xcas, odesolve supports solving first-order systems (https://www-fourier.ujf-grenoble.fr/~par...de592.html), and testing on my Prime confirms that it is implemented (although, as far as I can tell, undocumented).

To evaluate x(t1) given initial conditions x(t0)=A and x'(t0)=B, we would enter the following:

odesolve([x1, f(t,x0,x1)], t=t0..t1, [x0, x1], [A, B])

It will return a vector with two entries. The first is x0(t1) = x(t1), and the second is x1(t1) = x'(t1).


Example: (from https://sites.science.oregonstate.edu/ma...o_num.html): Find y(1) if y''=-y'+sin(t*y), y(0)=1, and y'(0)=2.

The differential equation reduces to the first-order system

y0' = y1
y1' = -y1 + sin(t*y0)

In the Prime, we enter

odesolve([y1, -y1+sin(t*y0)],t=0..1,[y0, y1], [1, 2])

Which returns the vector

[2.4504622693,1.18423606772]

This suggests y(1) = 2.4504622693 (and also y'(1)=1.18423606772). This is very far off from the answer of 4.1278 given by the source of the example, but Mathematica also gives 2.4504623, so I suspect this is right.

[parisse]

Exactly, you must rewrite your equation as a 1st order resolved system (in Cauchy-Lipschitz form).

[Wes Loewer]

(11-15-2023 03:27 AM)Brian Zilli Wrote:  I don't know if this qualifies as built-in since it does require some pre-processing, but this is one way to do it:

Suppose we have a second-order equation x''=f(t,x,x'). If we let x0=x and x1=x', we may transform this into a system of first-order equations:

x0' = x1
x1' = f(t,x0,x1)

Thanks. As an educational exercise I wrote a little program that uses that idea with some Euler variations and RK4. It made me wonder if there might be something already built-in, perhaps undocumented.

[Wes Loewer]

(11-15-2023 06:57 AM)parisse Wrote:  Exactly, you must rewrite your equation as a 1st order resolved system (in Cauchy-Lipschitz form).

Thanks.

Since desolve accepts 2nd Order Diff Eq's directly, I had wondered if it might be possible to do the same with odesolve.