Solve Differential Equations in Python

[@trevorsimpson8788]

Let me tell you something, I was looking for your videos panicking for time, but here you are explaining every little detail. That's what make a great tutorial. Thanks Prof.

[@apm]

I'm glad it was helpful. Additional tutorials are at apmonitor.com/pdc/index.php/Main/SolveDifferentialEquations

[@hengky753]

Best channel so far about python!!

[@apm]

Thanks!

[@frozenburrito9313]

My homework is done in 5 mins. Thank you!

[@ssrwarrior7978]

These videos are awesome. I modified your code to meet my needs and it helped a lot. I wish KZbin had a rating system. I would rate contents on your channel 10 out 10. Again, thank you.

[@emmanuelmasemola1014]

You have saved my research projects, thank you

[@AJ-et3vf]

Thank you so much for these lucid, thorough tutorials. Very helpful

[@ronaldoadrian6701]

Thanks a lot. I used to solve these functions with matlab It's really interesting to try Python. Easy to understand~

[@apm]

Here are more examples: apmonitor.com/pdc/index.php/Main/SolveDifferentialEquations and in Gekko: apmonitor.com/pdc/index.php/Main/PythonDifferentialEquations (even easier). I'm glad you enjoyed it.

[@mpja]

Another way to implement the step is to define a u(t) with arbitrary parameters and define these parameters in the args in odeint. Ploting the u(t) can be done by using the np.vectorize function

[@apm]

Great observation! Example #4 implements this approach as shown here apmonitor.com/pdc/index.php/Main/SolveDifferentialEquations

[@tylerriley]

What’s the point of the variable tspan @21:45. It was defined but never used. Was it supposed to replace t in odeint?

[@mariomuysensual]

i dont get it as well

[@GreenCrap24]

Very helpful, much better than the Simple/Complex video by you guys

[@apm]

Thanks for the comment and feedback. Here are additional tutorials that you may find useful: apmonitor.com/pdc/index.php/Main/SolveDifferentialEquations There are also gekko tutorials as an alternative to ODEINT: apmonitor.com/pdc/index.php/Main/PythonDifferentialEquations

[@kristianmr8629]

Hello Professor, I would like to thank you for the excellent content. I have one question regarding problem 4. Could we solve it as problem 3 ? I mean avoiding the 'for' loop and also insterting the 'u(t)' in an 'if' statement inside the function 'model'

[@apm]

Yes, that is possible. I introduced the 'for' loop method because we need this approach later in the course where the u(t) is calculated based on the y(t) value as a particular sample time. This is used in control such as PID control: apmonitor.com/pdc/index.php/Main/ProportionalIntegralDerivative

[@afsanehgharouni9201]

Hi, Thanks for the video and codes. In problem 4, when you want to get a step response, the input value of a given time is fixed and fed to model function. But odeint changes the step size and moves in time (back and forth, if you print(t) in model function) . For example, while it works on time span of [2.7, 2.8], it jumps also to t=6.3. During this process u[i] is fixed and it is zero even for t=6.3. Won't is possibly cause erroranous results? I'm not an expert in this area and details of how the solver works but this observation made me confused.

[@apm]

It isn't a problem because even though ODEINT is using t=6.3 to calculate the result, only the values at [2.7, 2.8] are reported as the solution. Any values calculated further out are discarded and only used to obtain the solution. There is an error tolerance for the solution so the solver continues working until the condition is satisfied so you can typically trust the solution.

[@musicarroll]

Very well done.

[@edholanda5897]

Professor how would I solve this system of first order edos numerically by plotting the graph for the different values ​​of (n). the derivatives are in relation to ha (r). a'/r = -e^2*v^2*(g^2 - 1) g' = - a*g/r given the boundary conditions a(0) = n a(inf)=0 g(0) = 0 g(inf)=1 where (e)=0.5 and (v)=1 are constant. please give a helping hand there, I looked for and did not find any problems like this on the matlab website. I'm from Brazil.

[@aksa101010]

For the last problem: Why didn't you use if/else for u function, as you've already done in the Problem 2? Thanks!

[@apm]

You could do it that way as well. The advantage of this way is that you can feed in many different inputs that would be too complex to implement as if statements.

[@ivanluthfiihwani9867]

Very interesting, Thank you Prof.

[@homazeinali5659]

Awesome. How could I solve d^2y/dx^2=y(x)/1+y(x) with boundary conditions y'(0)=0 and y(a)=1 in which a is a parameter?

[@apm]

Example 1b is similar to your problem: apmonitor.com/do/index.php/Main/DynamicOptimizationBenchmarks Here is information on second order systems: apmonitor.com/wiki/index.php/Apps/2ndOrderDifferential

[@saeedmardani3900]

Like always, you are amazing (y)

[@apm]

Thanks, Saeed!

[@ibrahimaba8966]

Hello. I am confuse. You define tspan but you didnt use it anywhere in the odeint function. I think you must replace t by tspan. Sorry for my english. From Paris.

[@apm]

Thanks, I fixed this in the code posted here: apmonitor.com/pdc/index.php/Main/SolveDifferentialEquations

[@ibrahimaba8966]

@@apm Thank you.... this video is interesting and helpful.

[@devrajsadaula4888]

Thank you Professor. Could you please give an example in solving system of linear differential equations where the parameter are also function of time. For example at time 28:30 of this video, kr1, kr2 .... are constant parameter. How do we solve if these are also function of time?

[@apm]

Example 4 shows how to change a parameter value during the integration: apmonitor.com/pdc/index.php/Main/SolveDifferentialEquations

[@CrushOfSiel]

Hello, I was wondering if any of your videos show code for graphing the direction field/isoclines? Thanks. I am really appreciating the content that I've seen so far. It's really helping me get back up to speed for this class I'm taking in graduate school.

[@apm]

I don't have anything but here is something on StackExchange that may help: stackoverflow.com/questions/18832763/drawing-directions-fields

[@CrushOfSiel]

Thank you! I'll check it out.

[@JackFrost-oo3dv]

Excuse me, can you suggest which book is the best for learning this one? Thank you

[@apm]

Online courses are available such as these: apm.byu.edu/prism/index.php/Site/OnlineCourses I don't recommend a book to learn programming because there are so many excellent online resources that are more interactive.

[@yem.t.3930]

Thank you very much! clear and concise!

[@lillianmulligan6521]

Are there units on time? can you specify weather it's minutes, hours, etc.?

[@apm]

Yes, you specify the units on the right-hand side of the equation. For example a mass balance is dm/dt = mflow_in - mflow_out. If the mflow_in and mflow_out are in kg/sec or kg/min then it will affect the solution for m(t).

[@universallappen1033]

Why does model take t as an argument if the return value of model is completely unrelated to t?

[@apm]

It tells how far to integrate in time. You only need the beginning and ending time at a minimum. If the model doesn't use time then t=[0,1] gives the same solution as t=[5,6].

[@quentinperret2981]

YOU ARE MY HERO, YOU JUST SAVE MY RESEARCH PROJECT! MARRY ME !!!

[@Iamduckpunk]

ok, may be a stupid question, but anyway what should i do if my condition are u(0) =0 and u'(1) = 1?

[@apm]

You can define an additional initial condition with ODEINT such as: import numpy as np from scipy.integrate import odeint import matplotlib.pyplot as plt # function that returns dy/dt def model(y,t): x,v = y dxdt = v # x = distance dvdt = 1.0 # v = velocity return [dxdt,dvdt] # initial condition y0 = [0,0] # time points t = np.linspace(0,5) # solve ODE y = odeint(model,y0,t) # plot results plt.plot(t,y) plt.xlabel('time') plt.ylabel('y(t)') plt.show() However, if you have a final condition, you'll need to use a package such as GEKKO: apmonitor.com/wiki/index.php/Main/GekkoPythonOptimization (see examples 11-13).

[@RakeshKumar-wx8hm]

Thanks for the nice lecture on python for differential equation. How to plot impulsive differential equations?

[@apm]

Numerically, a unit impulse is a step up and down with an integral of 1 (area). Ideally, the impulse goes to infinity for an infinitely small time step. Numerically, I suggest an impulse to a finite number such as 1,000,000 for 0.000001 time units.

[@MrAppleman666]

Cool video, thanks. Do you know how to force the odeint to deal with stiff problems? It seems that the automatic recognition of the problem's stiffness of the lsoda fortran code is not always working right.

[@apm]

Here are some suggestions: stackoverflow.com/questions/2088473/integrate-stiff-odes-with-python There is also the Gekko platform: gekko.readthedocs.io/en/latest/

[@virtualpeyman]

in problem 4 code why do you define tspan but do not use it in z definition in the loop?

[@apm]

I made a mistake when recording the video but it should give the same solution. Here is the corrected version: apmonitor.com/pdc/index.php/Main/SolveDifferentialEquations?action=sourceblock&num=6

[@virtualpeyman]

@@apm you are amazing!

[@komalkaursasan5425]

Pretty amazing. Where can I get the solution for the assignment question?

[@apm]

Here are the solutions: apmonitor.com/pdc/index.php/Main/SolveDifferentialEquations

[@komalkaursasan5425]

@@apm Thank you!

[@jdcrunchman999]

What does the variable “k” mean.. that was not explained

[@apm]

It is the rate constant that is related to the half-life if this were a radioactive decay problem.

[@penelopewaters4630]

I was looking for this comment thank you

[@csegura26]

Very helpful tutorial .... thanks!

[@ssrwarrior7978]

This is an Awesome video.. How to solve ODE if I have array of parameters that changes with time

[@apm]

Here is more help related to time varying inputs: apmonitor.com/pdc/index.php/Main/SolveDifferentialEquations

[@apm]

Here is another example on Stack Overflow: stackoverflow.com/questions/74362585/bioreactor-simulation-for-ethanol-production-using-gekko/74369410#74369410

[@ssrwarrior7978]

@@apm Thank you for the reply. I am trying to solve voltage response across RC circuit with three RCs in series.. I was trying to your method show in problem 4 but could not use v = odeint(v_response,v0,tspan,args=(r1[j],c1[j],r2[j],c2[j],r3[j],c3[j],i[j])) this gives error "v_response() missing 1 required positional argument: 'I' here 'I' is current array that changes with time, t. Also r1,c1; r2,c2; r3,c3; change with time ...Not sure why it his not working

[@ssrwarrior7978]

j is the loop counter

[@ssrwarrior7978]

This is the function: This is common ODE for solving battery voltage def v_response(v,v0,t,r1,c1,r2,c2,r3,c3,i): v1=v0[0] v2=v0[1] v3=v0[2] dv1dt=(i/c1)-(v1/(r1*c1)) dv2dt=(i/c2)-(v2/(r2*c2)) dv3dt=(i/c3)-(v3/(r3*c3)) dvdt=[dv1dt, dv2dt, dv3dt] return dvdt

[@camoph]

Hi. How can I plot the derivatives, for instance dx/dt, together with the solution x(t)?

[@apm]

The model function returns the derivatives. After you solve the problem, you can call the model function with values of the state (x) to report the derivative values (dx/dt).

[@dmkdata]

This was very helpful, thanks alot!

[@qsnipes91]

hello, I was wondering if you have any examples for second order ode for initial value problems?

[@apm]

Sure, here is some code for 2nd order systems: apmonitor.com/wiki/index.php/Apps/2ndOrderDifferential

[@apm]

Here is another link: apmonitor.com/pdc/index.php/Main/SecondOrderSystems

[@qsnipes91]

@@apm Thank you! I'll take a look

[@mpja]

I'm a bit confused at the function x(t) for example 4. If x(t) the integral of dxdt, why is it being treated as an independent function in the dydt equation? Shouldn't dy/dt = (-y(t) + integral (dx/dt))/5 then?

[@apm]

Example 4 at 19:50 is correct. The value x(t) is the solution to the differential equation, not the integral of dx/dt. The first equation could be used to solve for x and then substituted into the second equation to solve for y. However, ODEINT has no problems solving many simultaneous differential equations. You just need to give the initial conditions and a function that returns the derivatives. There is more information on ODE solutions here: apmonitor.com/che263/index.php/Main/PythonDynamicSim I recommend that you go through Euler's method as an introductory exercise to give insight on how an integrator (such as ODEINT) works. ODEINT uses higher order methods such as Runga-Kutta, however.

[@vijaykrishnamenon]

I tried this way and I get the same result, Why is it required to do integrate in multiple steps when this work. Did I do anything wrong??? This is just the same way the 2nd example was done. But I get it by just one call to ODEINT!! # solving 2 equations simultaneously ... # 1. 2.dx/dt = -x(t)+ u(t), x(0) = 0, # u = 2.S(t-5) and S(t-5) is a step # function which toggles to 1 when t=5 # 2. 5.dy/dt = - y(t)+ x(t), y(0) = 0 # call to the function that return dy/dt def model(z,t): if (t < 5 ): u = 0.0 else: u = 2.0 x = z[0] y = z[1] dxdt = (-x + u)/2.0 dydt = (-y + x)/5.0 return [dxdt,dydt] # initial condition z0 = [0,0] # time points n = 150 t = np.linspace(0,15,n) # solve ode z = odeint(model,z0,t) #plot resutls x= z[:,0] y= z[:,1] plt.plot(t,x,'r-',linewidth=2,label='x(t)') plt.plot(t,y,'b--',linewidth=2,label='y(t)') plt.xlabel('time') plt.ylabel('z(t)') plt.legend() plt.show()

[@apm]

This way is perfectly acceptable. I show how to break it into multiple steps for later exercises in the Process Dynamics and Control course where there is a feedback controller (PID) where the value of the input (u) is not known until you integrate forward one time step. This helps us use a process simulator to emulate discrete control such as apmonitor.com/pdc/index.php/Main/LevelControl

[@vijaykrishnamenon]

@@apm thank a lot. It is clear now. I need to solve for systems of first order ODEs... this too is very helpful as mine is data driven too.

[@apm]

Correction to the video at 22:00 - it should be tspan inside the loop instead of t in the call to odeint. It is correct on the web-site at apmonitor.com/pdc/index.php/Main/SolveDifferentialEquations (see Problem 4) and in the new video kzbin.info/www/bejne/nIPTYn9_rcyaetU

[@maryamasgharpour7910]

Here you ran tspan or x. I am confused. how can I run that? I use updated python in Mac not Anecdote.

[@apm]

@@maryamasgharpour7910 here is the code that you can copy: apmonitor.com/pdc/index.php/Main/SolveDifferentialEquations?action=sourceblock&num=6

[@niuer2]

all the engineer students should see this channel

[@apm]

Thanks wiaowen chen!

[@MultiGuitarman21]

When solving the ODE iteratively, in the ''for loop' can someone explain why "z0 = z[1]"? I mean shouldn't it be "z0=z[i]"

[@apm]

For every step, you take the last value to be your new initial condition. The result z is only length 2 so z[1] is the result from integrating from z0 (initial condition) to the next time step given by tspan. You could also write z0=z[-1] to always get the last value if tspan has 3+ values.

[@DipeshKumar-gc2mi]

@@apm z has length 2 then why only one boundary condition is updated i.e. z[1], why not both boundary condition

[@rrc]

@@DipeshKumar-gc2mi you would need to update z0=[z[-1,0],z[-1,1]] if you have two states.

[@sampaiomag1]

It is amazing. Thank you so much.

[@pnachtwey]

Lin space doesn’t feel right to me. I think in terms of number on intervals, not points. I am always adding one to the third parameter to make the time points nice numbers. I prefer changing u is the model so the odeint doesn’t need to be in a loop. There are incremental forms of PID so u is integrated too.

[@apm]

You may like the np.arange function to give you the intervals you desire.

[@natism224]

Good moring. Thanks for this video. What if I have two differential equations but one depends on the other one?

[@apm]

Natalia, please see problem #4 here: apmonitor.com/pdc/index.php/Main/SolveDifferentialEquations the equations are solved together. You can also see the Gekko tutorial link that covers those same problems. Gekko is a bit easier to understand than ODEINT.

[@muratkaradag3703]

File "", line 14 y == odeint(model, y0 , t) ^ SyntaxError: invalid syntax if I type that code ....

[@apm]

Just use one equal sign.

[@muratkaradag3703]

@@apm haha i found it, Kaggle showed the mistake in that line, but i didnt do a closing parentheses as i defined t :) THX FOR YOUR ANSWER

[@yaelio9215]

Thank you for this superb explanation. I wonder if you have an example for solving a set of ODEs where some of them contain more than one derivative in time, meaning something like: dx1/dt=… dx2/dt+c*dx1/dt=… could the second equation just be written as is into the function or some manipulations are required?

[@apm]

Please see Problems 3 and 4 for ODEINT (apmonitor.com/pdc/index.php/Main/SolveDifferentialEquations) or for GEKKO (apmonitor.com/pdc/index.php/Main/PythonDifferentialEquations) in Python.

[@yaelio9215]

Thank you very much but these examples unfortunately do not answer my question. Let's say that if my equations would be similar to those in Problem 3, only with slight difference in the 2nd equation: dxdt = 3.0 * np.exp(-t) dy/dt + dx/dt = -z[1] + 3 How can I tackle this? Because the second equation gives an error: "SyntaxError: can't assign to operator" so I wonder if there is a way to solve this in Python or I have to apply some mathematical manipulations to stay with one derivative in time for each equation.

[@rrc]

@@yaelio9215 You need to return the derivatives as a vector [dxdt, dydt]. You'll need to rearrange your second equation so that you only have the derivative of y on the left side. That way you can return the derivatives of X and Y from the function.

[@tonikpayc]

In 2:00, what would be the difference between making it model(y,t) and model(t,y)? Ive tried it and I get different numbers, so I know that there is, indeed, a difference. I would like to know what it is in specific.

[@apm]

The documentation for odeint shows that you need the states (y) listed before the time values (t). To see the documentation, run the following Python commands: from scipy.integrate import odeint help(odeint) MATLAB uses the other convention (t,y) so it can be confusing for those with experience with ode23 or other MATLAB integrators.

[@tonikpayc]

Thanks! That helped a lot. Just to avoid confusion though, the reason for this is simply because thats just the notation, correct?

[@apm]

Yes, that is just the way the developers require the order of the inputs.

[@zcq0730]

Very nice video, thanks!

[@hemantpandey9676]

How to fix "ModuleNotFoundError"..please help me !! I have recently installed python.

[@apm]

You may be missing SciPy. Here are instructions on installing packages: apmonitor.com/pds/index.php/Main/InstallPythonPackages such as SciPy and Numpy

[@hemantpandey9676]

@@apm Thanks..!

[@pr4wn5tar]

Top quality content. Thanks!

[@edsondepinhodasilva3696]

Great! Thank you very much!

[@moirangthemsanahal6695]

Thanks. Need another on 2nd orders

[@apm]

Sure, here is additional information on 2nd order systems: apmonitor.com/pdc/index.php/Main/SecondOrderSystems and apmonitor.com/wiki/index.php/Apps/2ndOrderDifferential

[@kasperlakshitha7768]

dx/dt=4(x^2+1), x(pi/4) = 1 Dude I'm stuck with this one. Help me out.

[@apm]

The first example here: apmonitor.com/pdc/index.php/Main/SolveDifferentialEquations should help. You can copy the code and put in your equation and initial condition.

[@vijaydhari6235]

Sir, I want to create a shape similar to the amoeba. The amoeba which I create should fit in between two concentric circles, the boundary condition what I am thinking is since in amoeba the starting point and ending point is the same and the shortest distance between two diametrically opposite points should be greater than 2r of the inner circle.

[@apm]

That sounds like an interesting application. I'm not sure why you requested that on this video but hopefully you can find some help online.

[@vijaydhari6235]

@@apm Sir I am working on reservoir simulation and any online resource would be of great help.

[@apm]

@@vijaydhari6235 Here is an additional resource that may help: apmonitor.com/pdc/index.php/Main/PythonDifferentialEquations

[@galalump9405]

Very informative.... thank you

[@apm]

Glad it was helpful!

[@kqnrqdtqqtttel1778]

PROBLEM 4 ------------------- def functionModel5(z2, t): if t < 5: u = 0 else: u = 1 x = z2[0] y = z2[1] dxdt = (-x + u)/2 dydt = (-y + x)/5 return dxdt, dydt z2_0 = [0, 0] t = np.linspace(0, 14) z2 = odeint(functionModel5, z2_0, t) x = z2[:, 0] y = z2[:, 1] plt.plot(t, x, 'b-') plt.plot(t, y, 'r:', linewidth = 2) plt.title("Step-function / double Function Solution", fontsize = 14) plt.xlabel('Time-Axis', fontsize = 13) plt.ylabel('Y-Axis', fontsize = 13) plt.legend(['x(t)', 'y(t)'], loc = 'best')

[@DanielLima-bd8sv]

Very helpful. Thanks. Could you please make videos on PDEs and sympy. That would be nice and of great help too.

[@apm]

Great suggestion. I haven't solved PDEs with Sympy. Here is a tutorial on solving PDEs with Python Gekko: apmonitor.com/do/index.php/Main/PartialDifferentialEquations Sympy is likely limited to very simple PDEs while numeric solutions, such as what Python Gekko produces, can handle much more complex problems.

[@nfkt101]

How did you obtain u[51:]=2 ?

[@apm]

At 21:08, the input (u) is adjusted to a value of 2 after the 5 seconds of simulation time. Because each step is 0.1, the first 5 seconds of time points are 0, 0.1, 0.2, ... 4.9, 5.0. The Python index for these values is 0,1,2, ..., 49, 50. We change the value of u after the first 5 seconds so starting at element 51. The syntax to change from 51 to the end to a value of 2 is u[51:] = 2. I hope this helps.

[@nfkt101]

@@apm Thanks!

[@penelopewaters4630]

what does the k mean?

[@apm]

It is the rate constant. Here are additional examples: apmonitor.com/pdc/index.php/Main/SolveDifferentialEquations

[@muhammadhamzashahid9649]

Thank you sir.

[@jadhavsourabh]

Thank you ☺️

[@AjeetKumar-oc8sb]

Hello sir, I'm getting error from scipy.intergate import odeint saying: ImportError: cannot import name 'ideint' Also I have scipy installed correctly

[@apm]

Make sure it is from scipy.integrate import odeint - you have some letters that aren't correct in integrate and odeint.

[@muhammadadeel1150]

Is this fine in case of ode45 ?

[@apm]

Yes, this is a good substitute for ode45.

[@muhammadadeel1150]

@@apm Thank you so much. May i have your email address?

[@mermichael9108]

How do I get the numeric solution?

[@apm]

Here are some examples on getting a numeric solution: apmonitor.com/pdc/index.php/Main/SolveDifferentialEquations If you want an analytic solution: github.com/APMonitor/data_science/blob/master/11.%20Differential_Equations.ipynb

[@mermichael9108]

@@apm I wanted to get the general solution for a differential equation. You can get that by using scipy package in Python. The method dsolve() returns the general solution.

[@apm]

@@mermichael9108 that is a good option for linear systems of differential equations. The ODEINT solver numerically solves both linear and nonlinear systems of differential equations. Another option is Python Gekko that solves Differential and Algebraic Equation (DAE) systems of higher index in simulation or optimization: apmonitor.com/pdc/index.php/Main/PythonDifferentialEquations

[@icabanah]

Ty so much!

[@franciscomiranda3239]

Grear man. Can you make a graphic voronoi hexagon report with x,y points datas??

[@apm]

Sorry, I don't know about that. Also, it doesn't seem related to this video. Maybe look on Stackoverflow for related content.

[@rahulbball9395]

how about boundary value problems using solve_bvp? any good tutorials for that?

[@apm]

Sorry, I don't have a tutorial on that function. Here is some documentation: docs.scipy.org/doc/scipy-0.18.1/reference/generated/scipy.integrate.solve_bvp.html

[@ameyakabra449]

include how to download the modules also please

[@apm]

apmonitor.com/pdc/index.php/Main/InstallPython

[@fgularte]

Python 3.6.4 numpy and first lines do not work :-L!

[@apm]

You may have a problem with the installation. Here are instructions on installing numpy, scipy, and matplotlib. Either method 1 or 2 should work. There are also other examples posted to apmonitor.com/pdc

[@apm]

Also, make sure you use the "get code" link in the bottom right of the source block or else you may get extra formatting characters. Here is the raw source. apmonitor.com/pdc/index.php/Main/SolveDifferentialEquations?action=sourceblock&num=1

[@xiuxiangwu2646]

@@apm I didn't help me at all. I just wanted to download python!

[@apm]

Here is the link with instructions on downloading and installing Python: apmonitor.com/pdc/index.php/Main/InstallPython This video is about solving differential equations.

[@jammutarun8900]

2nd order differential equation

[@apm]

Here is help on 2nd order systems: apmonitor.com/wiki/index.php/Apps/2ndOrderDifferential

[@godyt2548]

5:41 23:45

[@apm]

Did you have a particular question about this point in the video?

[@rickandelon9374]

Use Mathematica guys! Don't bother using python for the quality you get out Mathematica!

[@apm]

Mathematica is also a great program but a proprietary platform. One of the reasons that Python is so popular is for the full range of packages and capabilities. MATLAB is another proprietary package that I'd recommend over Mathematica.