Please I have the following problem: I have an equation for the elastoplastic response. I have split the beam into two segments (the first segment is defined by the coordinates (x1, y1) and (x2, y2) and the second segment is defined by (x1, y1) and (x2, y2).
Basically, the equations that come from the kinematics and the equilibrium forces and that I want to solve (numerically) (I am only interested in the situation i = 1, 2, 3)

$ dfrac {x_ {i + 1} + x_ {i}} {l} – cos theta_i = 0 $ (two equations to the system), $ theta_i $ is an angle in the ith beam connection.

$ dfrac {y_ {i + 1} + y_ {i}} {l} – sin theta_i = 0 $ (two equations to the system)

$ V_ {Xi} – V_ {Xi + 1} – F _ {Xi} = 0 $ (Forcing balance. $ F_ {Xi} $ denotes the external force and $ V_ {Xi} $ denotes the net force vector. It is assumed that there are no forces $ F $ or $ V $ in the direction of the xaxis.

$ V_ {Yi} – V_ {Yi + 1} – F _ {Yi} = 0 $ (One assumption is that there is no power $ F $ in the direction of the yaxis. $ V_ {Yi} $ are calculated so that it is simplified to the relationship $ V_ {Yi} = V_ {Yi + 1} $

$ M_ {i} – M_ {i + 1} + (V_ {Xi} – 1 / 2F_ {Xi}) * l * sin theta_ {i} – (V_ {Yi} – 1 / 2F_ {Yi}) * l * cos theta_ {i} = 0 $ (Balance of momentum – using the assumptions can be simplified to the form: 5a) $ M_ {i} – M_ {i + 1} – V_ {Yi} * l * cos theta_ {i} = 0 $ (two equations to the system). $ M_i $ is the bending moment at the beam connection.

Response relation:
$ dot {M_ {i}} = lbrace 1 – H left (M_i dot { kappa_i} right) left (1 + tanh left ( alpha H left (1 – left  dfrac {M_i} {M_y} right  right) left ( left  dfrac {M_i} {M_y} right  – 1 right) right) rbrace EI cdot dot { kappa_i} $ (An equation to the system – this answer relation holds only in the inner joints, in this case only i = 2). H denotes the Heaviside function, EI is the Young module, $ kappa $ is curvature, $ alpha $ is a coefficient, usually between 12 and 30 $ M_y $ is the flow moment that indicates whether the reaction is plastic or elastic (depending on whether $ M $ is less that $ M_y $).
I accepted the power $ Y_3 $ is charged by the stepbystep charging function $ f $, How $ f $ is not differentiable, I have its & # 39; derivative & # 39; manually defined as $ g $, So I have differentiated equation (5a) to include $ g $ (Hidden in everyone $ V_ {Yi} $).
For the Heaviside function $ H $ I used an approach $ H (x, n) approximately dfrac {1} {2} + dfrac {1} { pi} arctan (n * x) $, For more defined n, the approximation is sufficient.
My code is (including all initial and boundary conditions – the length of a segment $ l $ is 1 – I assume that only two segments of the total length are present $ 2 * l $, Also discretization for the curvature $ kappa $ how $ kappa_i approx dfrac { theta_i – theta_ {i1}} {l} $ is used.
n := 10000;
H(x_, n_) := (1/2) + (1/Pi) ArcTan(n*x);
My := 1; (*yield bending moment*)
(*approximation of the Heaviside step function*)
EI := 0.5; (*Young modullus*)
l := 1; (*length of one segment*)
tlimit := 1; (*for the definition of the loading function*)
Alpha1 := 30; (*coefficient*)
f(t_) := Piecewise({{t, t <= tlimit}, {tlimit,
t <= tlimit + 5}, {t  tlimit  6,
t <= tlimit + 6}}); (*loading function*)
g(t_) := Piecewise({{1, t <= tlimit}, {0, t <= tlimit + 5}, {1,
t <= tlimit +
6}}); (*derivative of loading function with definition of
problematic points*)
eq = {x2(t)  x1(t)  l*Cos(theta1(t)) == 0,
x3(t)  x2(t)  l*Cos(theta2(t)) == 0,
y2(t)  y1(t)  l*Sin(theta1(t)) == 0,
y3(t)  y2(t)  l*Sin(theta2(t)) ==
0, (theta2'(t)  theta1'(t))*(1 
H((theta2'(t)  theta1'(t))*M2(t)/l, n))*(1 +
Tanh(Alpha1*H(1  Abs(M2(t)/My), n)*(Abs(M2(t)/My)  1)))*(EI/
l) == M2'(t),
M1'(t)  M2'(t) ==
g(t)*l*Cos(theta1(t))  f(t)*l*Sin(theta1(t))*theta1'(t),
M2'(t)  M3'(t) ==
g(t)*l*Cos(theta2(t))  f(t)*l*Sin(theta2(t))*theta2'(t),
theta1(0) == 0,
theta2(0) == 0,
M1(0) == 0,
M2(0) == 0,
M3(0) == 0,
y1(0) == 0,
y2(0) == 0,
y3(0) == 0,
x1(0) == 0,
x2(0) == l,
x3(0) == 2*l,
x1(t) == 0,
y1(t) == 0,
M3(t) == 0,
theta1(t) == 0};
sol = NDSolve(eq, {x1, x2, x3, y1, y2, y3}, {t, 0, 100})
Plot(Evaluate(x3(t) /. sol), {t, 0, 20}, PlotRange > All)
Unfortunately, the code does not work. The error message is "Initial conditions that fulfill the residual function can not be found".
within the specified tolerances. Try to specify initial conditions for both.
Values and derivatives of the functions.
Can someone help me please ??? Thanks in advance!