2

I have a code for the solution of the Lane-Emden equation which is giving correct values after solving but when I plot those lists using matplotlib, I am getting blank images and not correct plots.

import numpy as np
import matplotlib.pyplot as plt

n = 14
theta_0 = 1
phi_0 = 0
h = 0.01
xi_0 = 0
xi_max = 100
   
theta = theta_0
phi = phi_0
xi = xi_0 + h
   
Theta = [[] for i in range(n)]
Phi = [[] for i in range(n)]
Xi = [[] for i in range(n)]
    
for i in range(n):
    Theta[i].append(theta)
    Phi[i].append(phi)
    Xi[i].append(xi)
     
def dThetadXi(phi,xi): #r1
    return -phi/xi**2
        
def r2(phi,xi):
    return dThetadXi(phi+h,xi+h*dThetadXi(phi,xi))
    
def r3(phi,xi):
    return dThetadXi(phi+h,xi+h*r2(phi,xi))
      
def r4(phi,xi):
    return dThetadXi(phi+h,xi+h*r3(phi,xi))
    
def dPhidXi(theta,xi,n): #k1
    return theta**(n)*(xi**2)
    
        
def k2(theta,xi,n):
    return dPhidXi(theta+h,xi+h*dPhidXi(theta,xi,n),n)
    
def k3(theta,xi,n):
    return dPhidXi(theta+h,xi+h*k2(theta,xi,n),n)
        
def k4(theta,xi,n):
    return dPhidXi(theta+h,xi+h*k3(theta,xi,n),n)
        
for i in range(n):  
    while xi < xi_max:
        if theta < 0:
            break
        dTheta = (step/6)*(dThetadXi(phi,xi)+2*r2(phi,xi)+2*r3(phi,xi)+r4(phi,xi))
        dPhi = (step/6)*(dPhidXi(theta,xi,i/2.)+2*k2(theta,xi,n)+2*k3(theta,xi,n)+k4(theta,xi,n))
        theta = theta+ dTheta
        phi = phi +dPhi
        xi = xi + h
        Theta[i].append(theta)
        Phi[i].append(phi)
        Xi[i].append(xi)
    print (i/2., round(xi,2), round(dThetadXi(phi,xi),2), round(xi/3./dThetadXi(phi,xi),2), round(1./(4*np.pi*(i/2.+1))/dThetadXi(phi,xi)**2,2))
    theta = theta_0
    phi = phi_0
    xi = xi_0 + h
    plt.plot(phi, xi)
    plt.show()

The correct plots should like this. Thanks.

EDIT After implementing the edits suggested below, I am getting this error output.

usr/local/lib/python3.7/dist-packages/numpy/core/_asarray.py:136: VisibleDeprecationWarning: Creating an ndarray from ragged nested sequences (which is a list-or-tuple of lists-or-tuples-or ndarrays with different lengths or shapes) is deprecated. If you meant to do this, you must specify 'dtype=object' when creating the ndarray
  return array(a, dtype, copy=False, order=order, subok=True)
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
TypeError: float() argument must be a string or a number, not 'list'

The above exception was the direct cause of the following exception:

ValueError                                Traceback (most recent call last)
<ipython-input-20-df5bae47b875> in <module>()
     62     phi = phi0
     63     xi = xi0 + step
---> 64     plt.plot(Phi, Xi)
     65 

7 frames
/usr/local/lib/python3.7/dist-packages/numpy/core/_asarray.py in asarray(a, dtype, order)
     81 
     82     """
---> 83     return array(a, dtype, copy=False, order=order)
     84 
     85 

ValueError: setting an array element with a sequence.

EDIT2 enter image description here

Isha Upadhyay
  • 23
  • 5
  • It seems like you're just plotting a single point here: `plt.plot(phi, xi)`. Shouldn't it be `plt.plot(Phi, Xi)`, with capital letters as the lists you initialized earlier? `plt.show()` should also be called after the for-loop. – Tinu Sep 06 '21 at 10:20
  • Oh yeah I see, my bad. But those corrections led to me having these errors which I have shown in the edits. – Isha Upadhyay Sep 06 '21 at 10:39
  • `TypeError: float() argument must be a string or a number, not 'list'` Your're tryting to plot a list of lists if I see that correctly. Convert it to `numpy` arrays or plot each sublist individually. – Tinu Sep 06 '21 at 10:44
  • I am not sure how to take a list out of the list of lists and plot it individually. (since I am new to Python.) – Isha Upadhyay Sep 06 '21 at 10:49
  • 1
    can I use '''plt.plot(Phi[0], Xi[0])''' as a plot of the first element? – Isha Upadhyay Sep 06 '21 at 10:50
  • For example `plt.plot(Xi[0],Phi[0])` if you like to plot the first sublist – Tinu Sep 06 '21 at 10:51
  • You`re welcome. I added an answer which summarizes our comments. Feel free to accept it. – Tinu Sep 06 '21 at 11:06
  • You are not implementing RK4. You are not applying the RK method for a coupled, vector-valued system. Your application of the slopes to the independent variable is straight wrong. // You made an error in the array construction and corrected it the wrong way, you should initialize `Phi=[phi0]` and then just `Phi.append(phi)` to get simple arrays of the values. – Lutz Lehmann Sep 06 '21 at 14:26
  • And in total the code appears to duplicate https://stackoverflow.com/questions/52330621/lane-emden-solutions-using-rk4-hard-coded/52342579#52342579 with all the same errors – Lutz Lehmann Sep 06 '21 at 14:29
  • 1
    I have greatly expanded the answer to the linked original code. – Lutz Lehmann Sep 16 '21 at 11:27

1 Answers1

2

To summarize the comments under the original post which led to the solution:

Consider calling plt.plot() and plt.show() only once after the for-loop and plot individual sublists of Xi and Phi with plt.plot(Xi[0],Phi[0]) instead of single points like you did here plt.plot(xi, phi).

Tinu
  • 2,432
  • 2
  • 8
  • 20
  • My doubt is solved regarding the original question, but after plotting the results the graphs doesn't resemble the ones linked in the question. I have implemented simple RK4 in the code and took the exact differetial equation to solve. I don't understand the problem. Would you be able to provide some feedback regarding that? – Isha Upadhyay Sep 06 '21 at 12:34
  • Unfortunately I have no knowledge about the underlying problem. I know how RK4 works, but I've never heard of the Lane-Emden equation. Chances are there's a bug in your RK4 implementation. Can you verify your implementation with a simpler ODE? – Tinu Sep 06 '21 at 12:45
  • Yeah sounds like a good idea. Thanks anyway for the pointers! – Isha Upadhyay Sep 06 '21 at 13:00