Here is the Lorenz Attractor for MATLAB. It runs in Octave but the plot step size is too large and the plots are not very neat. It works great in MATLAB. Fixing it for Octave is left as an exercise for the reader.
%Lorenz Equation
sig=10; b=8/3; r=24.5;
f=@(t,y)[sig*(y(2)-y(1)); r*y(1)-y(2)-y(1)*y(3);
y(1)*y(2)-b*y(3) ];
[t,y]=ode45(f,[0 50],[-50; 0; 50]);
%Plot of x and y versus t
figure(1)
plot(t,y(:,1),t,y(:,2))
title('9.3.7 Plot X and Y versus t')
xlabel('t')
ylabel('x and y')
legend('x','y')
%Plot of x and z
figure(2)
plot(y(:,1), y(:,3))
title('9.3.7 Plot of X versus Z')
xlabel('x')
ylabel('z')
figure(3)
plot3(y(:,1),y(:,2),y(:,3))
title('9.3.7 3D Plot')
xlabel('x')
ylabel('y')
zlabel('Z')
It is a system of 3 cross coupled ODEs:
https://en.wikipedia.org/wiki/Lorenz_system (https://en.wikipedia.org/wiki/Lorenz_system)
There is also an RLC circuit. Unzip the two files and run the RLCexample.m file. It will find the RLC.m function file.
For the Python crowd (I am using Spyder 4):
I have been working through a book "Deep Learning From Scratch ... " which gets deeply into partial differential equations. Here is the code for an example early in the book:
# -*- coding: utf-8 -*-
import matplotlib.pyplot as plt
import numpy as np
from numpy import ndarray
%matplotlib inline
from typing import Callable
from typing import List
Array_Function = Callable[[ndarray],ndarray]
Chain = List[Array_Function]
def square(x: ndarray) -> ndarray:
return np.power(x, 2)
def sigmoid(x:ndarray) -> ndarray:
return 1/(1+np.exp(-x))
def chain_length_2(chain: Chain,
x: ndarray) -> ndarray:
'''
Evaluates two functions in a row, in a "Chain".
'''
assert len(chain) == 2, \
"Length of input 'chain' should be 2"
f1 = chain[0]
f2 = chain[1]
return f2(f1(x))
def deriv(func: Callable[[ndarray],ndarray],
input_:ndarray,
delta: float = 0.001) -> ndarray:
return (func(input_ + delta) - func(input_ -delta))/(2*delta)
def chain_deriv_2(chain:Chain,input_range:ndarray)->ndarray:
assert len(chain) == 2,\
"This function requires 'Chain' objects of length 2"
assert input_range.ndim == 1,\
"Function requires a 1 dimensional ndarray as input"
f1 = chain[0]
f2 = chain[1]
f1_of_x = f1(input_range)
df1dx = deriv(f1,input_range)
df2du = deriv(f2,f1_of_x)
return df1dx * df2du
def chain_deriv_3(chain:Chain,input_range:ndarray)->ndarray:
assert len(chain) == 3, \
"This function requires 'Chain' objects of length 3"
f1 = chain[0]
f2 = chain[1]
f3 = chain[2]
f1_of_x = f1(input_range)
f2_of_x = f2(f1_of_x)
df3du = deriv(f3,f2_of_x)
df2du = deriv(f2,f1_of_x)
df1dx = deriv(f1,input_range)
return df1dx*df2du*df3du
def plot_chain(ax,
chain: Chain,
input_range: ndarray) -> None:
'''
Plots a chain function - a function made up of
multiple consecutive ndarray -> ndarray mappings -
Across the input_range
ax: matplotlib Subplot for plotting
'''
assert input_range.ndim == 1, \
"Function requires a 1 dimensional ndarray as input_range"
output_range = chain_length_2(chain, input_range)
ax.plot(input_range, output_range)
def plot_chain_deriv(ax,
chain: Chain,
input_range: ndarray) -> ndarray:
'''
Uses the chain rule to plot the derivative of a function consisting of
two nested functions.
ax: matplotlib Subplot for plotting
'''
output_range = chain_deriv_2(chain, input_range)
ax.plot(input_range, output_range)
fig, ax = plt.subplots(1, 2, sharey=False, figsize=(16, 8)) # 2 Rows, 1 Col
chain_1 = [square, sigmoid]
chain_2 = [sigmoid, square]
PLOT_RANGE = np.arange(-3, 3, 0.01)
plot_chain(ax[0], chain_1, PLOT_RANGE)
plot_chain_deriv(ax[0], chain_1, PLOT_RANGE)
ax[0].legend(["$f(x)$", "$\\frac{df}{dx}$"])
ax[0].set_title("Function and derivative for\n$f(x) = sigmoid(square(x))$")
ax[0].grid()
plot_chain(ax[1], chain_2, PLOT_RANGE)
plot_chain_deriv(ax[1], chain_2, PLOT_RANGE)
ax[1].legend(["$f(x)$", "$\\frac{df}{dx}$"])
ax[1].set_title("Function and derivative for\n$f(x) = square(sigmoid(x))$");
ax[1].grid()
It produces a couple of plots of chained functions and their derivative. Notice how the derivative is calculated. Just like the slope process from Calc I.
https://www.amazon.com/gp/product/1492041416 (https://www.amazon.com/gp/product/1492041416)
I omitted my preferred method of dealing with differential equations: MATLAB and SimuLink. Just plunk down integrators and gain blocks, wire them up and hit the "Go" button. This is the analog computer approach. Lord Kelvin rules!
Here's a MATLAB solution to a simple algebra problem (if solving 6 simultaneous equations is your cup of tea)
% A football stadium has a capacity of 100,000 attendees
% Ticket prices:
% Student = $25
% Alumni = $40
% Faculty = $60
% Public = $70
% Veterans = $32
% Guests = $ 0
% A full capacity game netted $4,987,000
% Let S = number of students, A = number of alumni, F = number of faculty
% P = number of public, V = number of veterans, G = number of guests
% S A F P V G = value
M = [ 1 1 1 1 1 1 100000; ... % total number of attendees
25 40 60 70 32 0 4897000; ... % total revenue
0 1 -1 0 0 0 11000; ... % 11000 more alumni than faculty
0 1 0 1 -10 0 0; ... % public + alumni = 10 * veterans
-1 1 1 0 0 0 0; ... % alumni + faculty = students
1 0 1 0 -4 -4 0]; % faculty + students = 4 (guests + veterans)
A = M(:,1:6); % extract 6x6 matrix
B = M(:,7); % extract value column vector
P = M(2,1:6); % ticket price by attendee type
R = inv(A)*B; % compute number of attendees by type
T = sum(R); % check total attendees = 100,000 CHECK
U = P*R; % total revenue = 4,897,000 CHECK
V = R'.*P; % revenue by attendee type - transpose R from 6x1 to 1x6
% then do element by element multiplication
%
% fancy output, all previous default output semicoloned
%
label = {'Students' 'Alumni' 'Faculty' 'Public' 'Veterans' 'Guests'};
for i = 1:6
fprintf('%-8s%7d @ %2d = %7d\n',string(label(i)),...
round(R(i)),...
M(2,i),...
round(V(i)))
end
fprintf('\nTotals %6d %8d\n',round(T), round(U))
Results:
Students 25000 @ 25 = 625000
Alumni 18000 @ 40 = 720000
Faculty 7000 @ 60 = 420000
Public 42000 @ 70 = 2940000
Veterans 6000 @ 32 = 192000
Guests 2000 @ 0 = 0
Totals 100000 4897000
A statement like
A = M(:,1:6); % extract 6x6 matrix
means fill a matrix A with all the rows of matrix M but only columns 1 through 6. I created M with 7 columns so I could document the output of each row. Then I break it into a 6x6 and 6x1 matrix to solve the equations.
It is much easier to start with something like this than to jump into differential equations. The syntax of MATLAB/Octave isn't all that simple. About the only way to learn it is to write it and slog through the errors.
There are hundreds of MATLAB books and I have a dozen or so. Pick one, any one, and work through some problems. Many books are of the "MATLAB For Physics" or other specific application.
Same story with Python - only far worse. How Python is taught is highly dependent on a supposed application. The Python syntax idea of using starting column to denote code level is wretched. Spaces shouldn't be a syntactic element. Nevertheless, Python is important, not because of the language, but because of the libraries.
Again, there are hundreds of books on Python and they too will tend to be application oriented. "PYTHON for Physics", that kind of thing. That "...Deep Learning..." book I linked above is an example.
As for coding stuff w/ Python vs a deicated math program, how hard is it to plot functions/etc w/ Python ? I suppose people have already made stuff for that.
Reply #9 has Python code (save as Fundamentals.ipy - the .i in .ipy is REQURED). This is a rather complex set of plots, not in the plotting, that is easy, but in the creation of the derivatives of chained functions.
Install Spyder 4 to use the code and after you press the dark green run button, the plots show up in the Plots tab over on the right of the screen.
Here is something from Google
Again, save this code as SineWave.ipy and run it inside Spyder 4
import matplotlib.pyplot as plt
import numpy as np
%matplotlib inline
Fs = 8000
f = 5
sample = 8000
x = np.arange(sample)
y = np.sin(2 * np.pi * f * x / Fs)
plt.plot(x, y)
plt.xlabel('sample(n)')
plt.ylabel('voltage(V)')
plt.show()
It doesn't get much easier than this! It probably works in every other IDE but I'm using just Spyder 4.
I went back and added the plot output from Fundamentals.ipy to Reply #9. It looks pretty nice considering how little effort goes into it.