Solving 2nd order differential equations wrt this code

Question

I cannot write the program which is solving 2nd order differential equation with respect to code I wrote for y'=y

I know that I should write a program which turn a 2nd order differential equation into two ordinary differential equations but I don!t know how can I do in Python.

P.S. : I have to use that code below. It's a homework

Please forgive my mistakes, it's my first question. Thanks in advance

from pylab import*
xd=[];y=[]
def F(x,y):
    return y
def rk4(x0,y0,h,N):
    xd.append(x0)
    yd.append(y0)
    for i in range (1,N+1) :
        k1=F(x0,y0)
        k2=F(x0+h/2,y0+h/2*k1)
        k3=F(x0+h/2,y0+h/2*k2)
        k4=F(x0+h,y0+h*k3)
        k=1/6*(k1+2*k2+2*k3+k4)
        y=y0+h*k
        x=x0+h
        yd.append(y)
        xd.append(x)
        y0=y
        x0=x
    return xd,yd
x0=0
y0=1
h=0.1
N=10
x,y=rk4(x0,y0,h,N)
print("x=",x)
print("y=",y)
plot(x,y)
show()

Already I ask how can I turn it into first order differential equation in python. Any example is enough for me, it’s not important what the differential equation is. By the way I canceled k4 as a mistake I will fix it. I’m sorry — user519955, May 28 '19 at 21:11
So you have `y''=f(x,y,y')` and rewrite it as `(y',v')=(v,f(x,y,v))`. Why should python turn it into a first order system? — Lutz Lehmann, May 28 '19 at 21:34
@LutzL OK. I’m trying to write code not to solve a 2nd order d.e. More clearly what should I add to my code for solving 2nd order d.e.? — user519955, May 28 '19 at 21:37
You need nothing to add, you need only make sure that all `y` values are always of `numpy` arrays containing the full state vector.. That is, `def F(x,y): return np.array([y[1], f(x,y[0],y[1])]);`. — Lutz Lehmann, May 28 '19 at 21:41
Related questions about using RK4 for second or higher order ODE: https://stackoverflow.com/q/52985027/3088138, https://stackoverflow.com/q/40919993/3088138, https://stackoverflow.com/q/53645649/3088138, https://math.stackexchange.com/q/1120984/115115, https://math.stackexchange.com/q/2615672/115115 — Lutz Lehmann, May 29 '19 at 07:09

score 0 · Answer 1 · answered May 28 '19 at 22:49

You can basically reformulate any scalar ODE (Ordinary Differential Equation) of order n in Cauchy form into an ODE of order 1. The only thing that you "pay" in this operation is that the second ODE's variables will be vectors instead of scalar functions.

Let me give you an example with an ODE of order 2. Suppose your ODE is: y'' = F(x,y, y'). Then you can replace it by [y, y']' = [y', F(x,y,y')], where the derivative of a vector has to be understood component-wise.

Let's take back your code and instead of using Runge-Kutta of order 4 as an approximate solution of your ODE, we will apply a simple Euler scheme.

from pylab import*
import matplotlib.pyplot as plt

# we are approximating the solution of y' = f(x,y) for x in [x_0, x_1] satisfying the Cauchy condition y(x_0) = y0

def f(x, y0):
    return y0

# here f defines the equation y' = y

def explicit_euler(x0, x1, y0, N,):
    # The following formula relates h and N
    h = (x1 - x0)/(N+1)

    xd = list()
    yd = list()

    xd.append(x0)
    yd.append(y0)

    for i in range (1,N+1) :
        # We use the explicite Euler scheme y_{i+1} = y_i + h * f(x_i, y_i)
        y = yd[-1] + h * f(xd[-1], yd[-1])
        # you can replace the above scheme by any other (R-K 4 for example !)
        x = xd[-1] + h

        yd.append(y)
        xd.append(x)

    return xd, yd

N = 250
x1 = 5
x0 = 0
y0 = 1
# the only function which satisfies y(0) = 1 and y'=y is y(x)=exp(x).
xd, yd =explicit_euler(x0, x1, y0, N)
plt.plot(xd,yd)
plt.show()
# this plot has the right shape !

Note that you can replace the Euler scheme by R-K 4 which has better stability and convergence properties.

Now, suppose that you want to solve a second order ODE, let's say for example: y'' = -y with initial conditions y(0) = 1 and y'(0) = 0. Then you have to transform your scalar function y into a vector of size 2 as explained above and in the comments in code below.

from pylab import*
import matplotlib.pyplot as plt
import numpy as np

# we are approximating the solution of y'' = f(x,y,y') for x in [x_0, x_1] satisfying the Cauchy condition of order 2:
# y(x_0) = y0 and y'(x_0) = y1

def f(x, y_d_0, y_d_1):
    return -y_d_0

# here f defines the equation y'' = -y

def explicit_euler(x0, x1, y0, y1, N,):
    # The following formula relates h and N
    h = (x1 - x0)/(N+1)

    xd = list()
    yd = list()

    xd.append(x0)
    # to allow group operations in R^2, we use the numpy library
    yd.append(np.array([y0, y1]))

    for i in range (1,N+1) :
        # We use the explicite Euler scheme y_{i+1} = y_i + h * f(x_i, y_i)
        # remember that now, yd is a list of vectors

        # the equivalent order 1 equation is [y, y']' = [y', f(x,y,y')]
        y = yd[-1] + h * np.array([yd[-1][1], f(xd[-1], yd[-1][0], yd[-1][1])])  # vector of dimension 2
        print(y)
        # you can replace the above scheme by any other (R-K 4 for example !)

        x = xd[-1] + h  # vector of dimension 1

        yd.append(y)
        xd.append(x)

    return xd, yd

x0 = 0
x1 = 30
y0 = 1
y1 = 0

# the only function satisfying y(0) = 1, y'(0) = 0 and y'' = -y is y(x) = cos(x)
N = 5000
xd, yd =explicit_euler(x0, x1, y0, y1, N)
# I only want the first variable of yd
yd_1 = list(map(lambda y: y[0], yd))

plt.plot(xd,yd_1)
plt.show()

Solving 2nd order differential equations wrt this code

1 Answers1