I have four ordinary differential equation written in mathmatica :

Simplify[DSolve[{y1'[x] == -0.162y1[x], y2'[x]=-0.148 y2[x]+ 0.055y1[x],y3'[x]==-0.134 y3[x]+ 0.033y1[x]+ 0.039y2[x],y4'[x]==-0.125y4[x]+0.021 y1[x]+0.025y2[x]+ 0.043y3[x], y1[0] == 100,y2[0]==4.76,y3[0]==69.7,y4[0]==0}, {y1[x], y2[x],y3[x],y4[x]}, x]]

I am trying to write them in python I have try the following :

def n1(y,x):
    s1 = 0.162131
    dydx = -s1 * y
    return dydx

def n2(y,x):
    s2 = 0.148234
    s12= 0.0555877 
    dydx2 = -s2 * y + s12 * n1(y,x)
    return dydx2

def n3(y,x):
    s3 = 0.134337
    s13= 0.0333526 
    s23= 0.0389114 
    dydx3 = -s3 * y + s13 * n1(y,x) + s23 * n2(y,x) 
    return dydx3

def n4(y,x):
    s4= 0.125072
    s14= 0.0213086
    s24= 0.0254777
    s34= 0.0430805
    dydt4 = -s4 * y + s14 * n1(y,x) + s24 * n2(y,x)  + s34 * n3(y,x) 
    return dydt4

## initial conditions
y01 = 100.
y02 = 4.74
y03 = 69.4
y04 = 0.    
x = np.linspace(0,20)

sol = odeint(n4,[y01,y02,y03,y04],x)

I want to solve n4 which include n1,n2,n3 , But I am only getting a solution for single equation??? I do not know what is wrong/missing to get the wright solution
Any help I appreciate that!!!

this calculation I can easly do it by hand. But I could not find a way to do it in python ?
I appreciate your help

I am trying to calculate some thing called directional derivatives of a vector b (which defined in attached file as an example)
to calculate the directional derivatives of the vector b .
the result of the example is also in the attached file.

I serached in goolge to find how we can calculate the directional derivatives, but I could not find any function in python to do this kind of calculation ; the only thing I found which do derivative in python is:

scipy.misc.derivative(func, x0, dx=1.0, n=1, args=(), order=3)

this function dose not give me the right answer .

Could you please help me to solve this probelm.
the example

this what I have tried to do , but I am still missing somthing which I do not know what is it?

def weighted(x, cols, w="weights"):
             return pd.Series(np.average(x[cols], weights=x[w], axis=0), cols)

 g.apply(weighted, ["column1", "column2"])

How can I calculate the average weight in python for two columns (x,y) data file, between 8<r< 9. in python

I have tried this put it dose not work with my

np.average(X,Y, weights=[8,9])

I have calculated the average weight which is 12.14, now I am trying to use this weight to calculate the average fro the X, Y values??
How can I do so ?

ant hints,sugesstion please ?

I have two array data files. I want to add weight function to these two array which I ploted see attached file. what I am trying to do now is use the weight function below, for two data set and select the points between 8<r <9 in each data then plot the histogram. how can I do so.

# weight function:
def W(r):
    Rsn = 5.4 
    return np.exp(-r/Rsn)

thank you very much

I solved it

plt.figure()
plt.hist(x[7]+ 1 ,bins=800,color='y',histtype='step',linewidth=2,stacked=True)
plt.show()

thanks for your help.

Any help please how to shift the plot to the right in the x-axis? I have tried put I did not find the right way to do so?

I want to shift my plot to the right by fixed value , I have npy file contain 7 columns, I plot the histogram for the column number 7 as below:

plt.figure()
plt.hist(x[7],bins=800,color='y',histtype='step',linewidth=2,stacked=True)
plt.show()

I am trying to shift my plot to the right with fix value but I could not find a correct way to do so. Any help I really appreciate it.

Thank you @xjr.If I want to add one parameter to the function and save the out put(sec_nucli) in a text file for example, I have try this:

Sec_nucli = 0
def Fragmentation(v,D,xsec):
    sec_nucli = v * n_H * D * xsec
    return sec_nucli

file=open("Secondary.txt", "a")     
for xi in xsec: 
    for Dj in D:
       print(xi, Dj)
     for vj in v:
       Sec_nucli += Fragmentation(vi,Dj, xi)
       print(Sec_nucli)
       s = str(Fragmentation(xi))
       file.write(s + "\n")
file.close()

Dose this correct!! .
Thank you very much for your help

I have tried this :

x     = np.array([16.7, 16.5, 16.1, 15, 13.5,14.2, 14.9, 13.7])# an example
D      = np.array([0, 1, 4, 7, 16, 31, 64, 127])# example
n_H = 0.9
Sec_nucli = []
x         = []
D         = []

for i in range(1,n): 
      Sec_nucli[i] =  n_H * D * xsec[i]
      sec_nucli += sec_nucli[i]
      return sec_nucli
## I defined the function like this:
def Fragmentation():
    n_H = 0.9
    Sec_nucli[i]=  n_H *D* x[i]
    return Sec_nucli

I want to define a function like this:

n_H  = 0.9 # constant . 
x     = np.array([16.7, 16.5, 16.1, 15, 13.5,14.2, 14.9, 13.7])# example to test
D      = np.array([0, 1, 4, 7, 16, 31, 64, 127])

Sec_nucli = []
x         = []
D         = []
Sec_nucli[i]=  n_H *D[i]* x[i]

where x, D are variables come from another functions, I put it as an array(just as a test here!!).
What I want is to make a for loop for the function Sec_nucli, that every time I use deferent value for x,D, and Sec_nucli. For that reson I put [i], infront of these variable(I dont know how can I write it)?!!. Could you please help me to write the function correctly?
Thank you for your help.

on the other word I have

  x=[x1,x2,....]
    x1: is the first point =
    ([27.266378342786531], [1.9505463431334713]),
    and so on..

how can I defined x ?

I have this point in (x,y) :
([27.266378342786531], [1.9505463431334713])
Now I want to use this point(I have many points,I want to know how can I use this point first ) to put it in the function below insted of x :(I have tried to make an array put i don't know how can I write it ??

def fun(x0, x, R, R0):
    x0=3.10e-24
    R= 3.15
    R0= 7
    return x0*np.exp(-np.maximum(x-R, 0)/R0)

x = np.array(???,???)
y= fun(x0, r, R, R0)
plt.plot(x, y)
plt.show()

thank you for your help

Thank you very much

I am trying to plot the picewise function below:

  rho0 = 3*10e-11  
  Rd = 3.15
  Rc = 7

    def Density(rho0, r, Rc, Rd):        
        if r < Rc :
          return  rho0 

        elif r> Rc :
          return rho0*np.exp(-(r-Rc)/Rd)

   r = np.linspace(0,20, 1000)
   y= Density(rho0, r, Rc, Rd)
   plt.plot(r, y)
   plt.show()

I am getting this error:
The truth value of an array with more than one element is ambiguous. Use a.any() or a.all().
I do not know how can I solve this problem?or

Thank you very much

I want to plot this function

H0*np.exp(r/Rh)
I definded it, then I want to plot it, r is different values ? yes you are right I will get exponential curve?

Okay, what about the plot how can I defined r

Hello,
I'm trying to plot this function with using this values:

Rh = 9.8    
H0 = 0.063 

def Thickness_H(H0, r, Rh): 
      return float(H0)*np.exp(float(r)/float(Rh))

r = np.linspace(0,1, 10)    
x =Thickness_H(H0, r, Rh)
plt.plot(r, x)
plt.show()

I'm trying to defind r like what i know but this dose not work, I'm getting this error:
only length-1 arrays can be converted to Python scalars.

may you please help me in this point?
Thank you

you mean like this:

float distance(float x, float y, float pitch)
int main()
{
float r,x,y,Bangle,pitch;  
      r = sqrt(x * x + y * y);
      Bangle = atan2(pitch*(y/r) + x , pitch*(x/r) -y);
      return (cos(Bangle),sin(Bangle),pitch);

}

I am trying to define a function in c++,
I wrote it in main function.
I am getting an error : a function-definition is not allowed here before ‘{’ token}I have defined it as following :

float distance(float x, float y, float pitch)
     {
      float r;  //local variable
      r = sqrt(x * x + y * y);
      Bangle = arctan2(pitch*(y/r) + x , pitch*(x/r) -y);
      return (cos(Bangle),sin(Bangle));
     }

I dont know what dose this means. Any sugestion

I am trying to define the archimedean spiral: when im trying to define the inclination angle (incl) of the tangent vector to the orbit ( i.e: tan(incl)) im getting an error :"'numpy.ufunc' object does not support item assignment".
the same error when I want to calculate cos(incl), and sin(incl).
Any suggestion and help.
My code is:

a= 2.     # will turn the spiral
v= 0.23
omega = 0.2
r0 = v/omega
t=np.linspace(0,T,N+1)
r= v*t
#theta = a + r/r0
theta  =  omega*t

x=r* np.cos( omega*t) 
y=r* np.sin( omega*t) 

dxdr = np.cos(theta) - (r/r0)*np.sin(theta)
dydr = np.sin(theta) + (r/r0)*np.cos(theta)

dydx = (r0*np.sin(theta) + r*np.cos(theta))/r0*np.cos(theta) - r*np.sin(theta)

np.tan[incl]= dydx
incl = np.arctan((dydx)) 

### Calculate cos(incl) ,sin(incl) :
np.sin[np.incl] = np.tan(np.incl)/np.sqrt(1+ np.tan(np.incl)*2)
np.cos[incl] = 1/np.sqrt(1+ np.tan(incl)*2)

Thank you for your explanation.
I have used generates Pseudo random numbers like this :

 Random random(0,1);   
    random.seed(2);   

It is work with me,thanks anyway

Well done.

or like this:

std::srand(std::time(NULL)); 
int random1 = std::rand();
int random2 = std::rand();

which one you think is better/correct

Thank you for your reply , I am using c++
I wrote it like this :

srand(time(NULL)); 
float   Random1 = (float(rand())/(RAND_MAX)) + 1
float   Random2 = (float(rand())/(RAND_MAX)) + 1
x = a* random1 + b* random2 
y = c* random1 + d* random2

is this correct ?

How can we insert two independent random number to two equations (x,y) contain two parts like this:
x = a randomnumber1 + b randomnumber2
y = c randomnumber1 + d randomnumber2
these two random number should be independent.
Thank inadvance.

I plot the archimedaen spirlas for one branch, how can I plot more than one branches starting from the original, my code that I have written :

n= 1000
a= 2.     # will turn the spiral
v= 1
omega = 0.1
r0 = v/omega
T=100
N=10000
t=np.linspace(0,T,N+1)
r= v*t
theta = a + r/r0
dt=float(T)/N 
x=r* np.cos(theta) 
y=r* np.sin(theta)
plt.plot(x,y)
plt.show()