If you are doing this for a coding product, a good possibility is doing it by numerically integrating the applied load $w$. The advantage of this is that it can be easily extended to general forms of $w$ for which closed expressions may be hard to get, the disadvantage is that it's more computationally expensive than evaluating a known formula.
The relevant equations are as follows, for the simply supported case shown in the image.
Reactions
$$
\begin{aligned}
P_2&=\frac{1}{L}\int_0^Lw(x)x\,dx \\
P_1&=\int_0^Lw(x)\,dx -P_2
\end{aligned}
$$
Internal shear and bending moment
$$
\begin{aligned}
V(x)&=\int_0^xw(x)\,dx-P_1 \\
M(x)&=-\int_0^xV(x)\,dx
\end{aligned}
$$
And, assuming Euler-Bernoulli beam theory,
$$
\begin{aligned}
\phi(x)&=\frac{M(x)}{EI}\\
\theta(x)&=\int\phi(x)\,dx +C \\
\delta(x)&=\int\theta(x)\,dx +Cx+D
\end{aligned}
$$
For the simply supported case $D=0$ and $C=-\left[\int\theta(x)\,dx\right](L)$
A MATLAB implementation of this is shown below
clear variables
close all
clc
%Inputs
%--Beam properties--
L=10; %Length, meters
E=200; %Modulus of elasticity, GPa
I=1; %Moment of inertia about axis considered, m^4
%--Load properties--
a=3; %Start of triangular load, m
b=4; %Length of triangular load,m
w=1; %Maximum magnitude, kN/m
%Boundary conditions (Only isostatic cases considered in this code)
fixed_slope=[0 0]; %0 indicates free dof, 1 means fixed dof.
fixed_deflection=[1 1];
%Simulation parameters
delta=0.01; %Max element length, m
%Calculations
% ADD A BIT HERE CHECKING THAT fixed_slope+fixed_deflection = 2 AND
% fixed_deflection has at leas one element or the problem is ill-posed
% ADD A BIT HERE CHECKING THAT a+b<=L and other validations
%
n=ceil(L/delta);
x=linspace(0,L,n)';
delta=x(2)-x(1); %Recalculated because original delta may have changed
q=(x>a).*(x<=a+b).*(w/b*(x-a)); %Conditional definition of applied load
% Equivalently using for loops and conditionals
% q=zeros(length(x),1);
% for i=1:length(x)
% if x(i)>=a && x(i)<=a+b
% q(i)=w/b*(x(i)-a);
% end
% end
figure
plot(x,q); xlabel('x [m]'); ylabel('Applied load [kN/m]')
%Determine reactions
qV=delta*trapz(q); %Total force applied by w
qM=delta*trapz(q.*x); %Moment generated by w around the leftmost end
P1=0; %Leftmost reaction, kN
P2=0; %Rightmost reaction, kN
M1=0; %Leftmost moment, kN-m
M2=0; %Rightmost moment, kN-m
% Five possibilities for reactions
if fixed_slope(1)==1 && fixed_deflection(1)==1
P1=qV; M1=qM;
end
if fixed_slope(1)==1 && fixed_deflection(2)==1
P2=qV; M1=qM-qV*L;
end
if fixed_slope(2)==1 && fixed_deflection(1)==1
P1=qV; M2=qM;
end
if fixed_slope(2)==1 && fixed_deflection(2)==1
P2=qV; M2=qM-qV/L;
end
if fixed_deflection(1)==1 && fixed_deflection(2)==1
P2=qM/L; P1=qV-P2;
end
%Calculate internal shear and bending moment
V=zeros(length(x),1); %Shear, kN
M=zeros(length(x),1); %Bending moment, kN-m
for i=1:length(x)
V(i)=delta*trapz(q(1:i))-P1;
end
figure
plot(x,V); xlabel('x [m]'); ylabel('Shear [kN]')
for i=1:length(x)
M(i)=-delta*trapz(V(1:i))-M1;
end
figure
plot(x,-M) %I plot moments on tensioned fibre side, just a convention
xlabel('x [m]'); ylabel('Bending moment [kN-m]')
%Some unit transformation, due to the way inputs were stated.
E=E*1e6; %kPa, or kN/m2
%Euler-Bernoulli beam theory calculation
phi=M./(E.*I);
theta=zeros(length(phi),1);
def=zeros(length(phi),1);
%Integrate radius of curvature to get slope
for i=1:length(theta)
theta(i)=delta*trapz(phi(1:i));
end
%Integrate slope to get deflection
for i=1:length(def)
def(i)=delta*trapz(theta(1:i));
end
% Five cases for boundary conditions
if fixed_slope(1)==1 && fixed_deflection(1)==1
C=-theta(1); D=-def(1);
end
if fixed_slope(1)==1 && fixed_deflection(2)==1
C=-theta(1); D=-def(end)-C*L;
end
if fixed_slope(2)==1 && fixed_deflection(1)==1
C=-theta(end); D=-def(1);
end
if fixed_slope(2)==1 && fixed_deflection(2)==1
C=-theta(end); D=-def(end)-C*L;
end
if fixed_deflection(1)==1 && fixed_deflection(2)==1
D=-def(1); C=-def(end)/L;
end
%Recall that we integrated twice, so we have a constant of integration and
%a constant of integration multiplied by x.
def=def+C*x+D;
figure
plot(x,def); xlabel('x [m]'); ylabel('deflection [m]')
%Plot a black line showing the original shape
hold on
plot([0 L],[0 0],'k','linewidth',2)
This produces the following results:
Applied load

Internal shear

Bending moment

Deflection

If you have any questions, just leave a comment and I'll get to it.