CODEWITHSUNDEEP
pythonintegration
Anne Sophie
Anne Sophie

Reputation: 187

Is there a way to integrate and get an array or a function instead of all the area under the curve?

I want to integrate the following equation:

d^2[Ψ(z)] / dz^2 = A * ρ(z)

Where Ψ (unknown) and ρ (known) are 1-D arrays and A is a constant.

I have already performed a Taylor expansion, i.e.

d^2[Ψ(z_0)] / dz^2 = [Ψ(z0+Δz) - 2Ψ(z0) + Ψ(z0-Δz)] / [Δz^2]

And successfully solve it by building a matrice.

Now, I would like to know if there is a Python (preferably) or Matlab function that can solve this function without having to do a Taylor expansion.

I have tried numpy.trapz and scipy.integrate.quad, but it seems that these functions only return the area under the curve, i.e. a number, and I am interested to get an array or a function (solving for Ψ).

Upvotes: 0

Views: 551

Answers (2)

Diego Palacios
Diego Palacios

Reputation: 1144

What you want to do is to solve the differential equation. Because it is a second order differential equation, you should modify your function to make it a system of first order ODEs. So you have to create a function like this:

Assuming ρ=rho

def f(z, y):
   return np.array([y[1], A*rho(z)])

where y is a vector containing Ψ in the first position and its derivative in the second position. Then, f returns a vector containing the first and second derivatives of Ψ.

Once done that, you can use scipy.integrate.solve_ivp to solve the problem:

scipy.integrate.solve_ivp(f, [z_start, z_final], y0, method='RK45', z_eval)

where y0 are the initial conditions of Ψ (the value of Ψ and its derivative at z_start). z_eval is the points where you want to store the solution. The solution will be an array containing the values of Ψ and its derivative.

Upvotes: 1

jeremy_rutman
jeremy_rutman

Reputation: 5788

You could integrate rho twice using an indefinite integral . e.g. 

import numpy as np
x=np.arange(0,3.1,0.1)
rho = x**3  # replace your rho here....

indef_integral1 = [np.trapz(rho[0:i],x[0:i]) for i in range(2,len(rho))]
indef_integral2 = [np.trapz(indef_integral1[0:i],x[0:i]) for i in range(2,len(indef_integral1))]

Upvotes: 0

Related Questions