Here is a Python implementation to generate shear and bending moment diagrams for statically determinate beams. It takes an arbitrary number of concentrated forces and the location of two reaction forces and then calculates the reaction force values and plots the shear and bending moment diagrams.

beam_1_diagram-a_darkversion-b

script_figure_output

import numpy as np
import operator
import matplotlib.pyplot as plt


### Input Parameters
#  l         Beam width/length.
#  r_R       Positions of supports (requires 2).
#  F_c       List of concentrated forces.
#              Enter in as a tuple the location measured from the leftmost part of the
#              beam, and the force value as the second item in the tuple. Upward is positive.

l =   30
r_R = [ 2, 15 ]
F_c = [ (5,-200), (13,-357), (20,75), (27,-147) ]


## Calculate reaction forces
moment_prod, force_sum = [], 0
for i in range(len(F_c)):
    moment_prod.append(F_c[i][0]*F_c[i][1])
    force_sum = force_sum + F_c[i][1]
moment_sum = sum(moment_prod)

A = np.array([[1,1],[r_R[0],r_R[1]]])
b = np.array([[-force_sum],[-moment_sum]])

R = np.linalg.solve(A,b)
R1 = float(R[0])
R2 = float(R[1])


## Do calcs for the shear diagram
# First create a load array and sort by location from left to right
loads_all = [ (r_R[0],R1), (r_R[1],R2) ]
for each in F_c:
    loads_all.append(each)

pos_all, f_all = [], []
loads_srtd = sorted(loads_all, key = operator.itemgetter(0))
for each in loads_srtd:
    pos_all.append(each[0])
    f_all.append(each[1])


# Calculate shear point values
V_pts = []
for i in range(len(f_all)):
    V_pts.append(sum(f_all[:i]))


# Create a shear value function
V_x, V = [], []
for i in range(0,len(V_pts)):
    V.append(V_pts[i])
    V.append(V_pts[i])
    V_x.append(pos_all[i])
    V_x.append(pos_all[i])

V.remove(0)
V.append(0)


## Do calcs for the moment diagram
# Calculate moment points
M_pts, M = [], []
for i in range(len(V)):
    M_pts.append(np.trapz(V[i:i+2],V_x[i:i+2]))

while 0 in M_pts:
    M_pts.remove(0)

for i in range(len(M_pts)):
    M.append(sum(M_pts[:i]))


# Create a moment function
M.append(0)
M_x = pos_all


## Print results
print(''.ljust(20),'Left'.ljust(20),'Right'.ljust(20))
print('Reaction Forces'.ljust(20),str(R1).ljust(20),str(R2).ljust(20))


## Create plots
# Plot a shear diagram
fl = 0.05*l
plt.subplot(211)
plt.plot([0,l],[0,0],'k-',linewidth=3)
plt.plot(V_x,V)
plt.xlim(-fl,l+fl)
plt.title('Shear Diagram')
plt.xlabel('x')
plt.ylabel('Shear Force')

# Plot a moment diagram
plt.subplot(212)
plt.plot([0,l],[0,0],'k-',linewidth=3)
plt.plot(M_x,M)
plt.xlim(-fl,l+fl)
plt.title('Moment Diagram')
plt.xlabel('x')
plt.ylabel('Moment')
plt.show()

Output window:

 
                     Left                 Right               
Reaction Forces      101.92307692307668   527.0769230769231