Invalid initial condition error

I have to solve the sistem of differential equation odesys with the condition imposed in bc vector. I obtain the "Invalid Initial Condition" at the line where v is defined, even if the domain for the boundary condition is correct. I must keep it a symbolic solution and a0 is a costant.

%% ANALYTICAL MODEL FOR A DCB SPECIMEN UNDER THE CONDITION OF PRESCRIBED DISPLACEMENTS
%% Linear Elastic Phase
%---------
syms x d v0(x) v1(x) v2(x) Lcz 
%---------
phi0 = -diff(v0,x);
M0 = E*I*diff(v0,x,2);
T0 = E*I*diff(v0,x,3);
phi1 = -diff(v1,x);
M1 = E*I*diff(v1,x,2);
T1 = E*I*diff(v1,x,3);
phi2 = -diff(v2,x);
M2 = E*I*diff(v2,x,2);
T2 = E*I*diff(v2,x,3);
%---------
ode_0 = diff(v0,x,4) == 0;
ode_1 = diff(v1,x,4) - 2*w*(lambda^2)*diff(v1,x,2) + (lambda^4)*v1 == 0;
ode_2 = diff(v2,x,4) + 2*ps*(k^2)*diff(v2,x,2) - k^4*(v2 - d_c/2) == 0;
%---------
syms xL xR xI
xL = -a0 - Lcz;
xI = -Lcz;
xR = L - a0 - Lcz;
c1 = v0(xL) == d/2;
c2 = M0(xL) == 0;
c3 = v0(xI) == v2(xI);
c4 = phi0(xI) == phi2(xI);
c5 = M0(xI) == M2(xI);
c6 = T0(xI) == T2(xI);
c7 = v1(0) == v2(0);
c8 = phi1(0) == phi2(0);
c9 = M1(0) == M2(0);
c10 = T1(0) == T2(0);
c11 = v1(xR) == 0;
c12 = phi1(xR) == 0;
%---------
odesys = [ode_0; ode_1; ode_2];
bc = [c1; c2; c3; c4; c5; c6; c7; c8; c9; c10; c11; c12];
v = dsolve(odesys, bc);
%---------
v1_sol(x,d,Lcz) = simplify(v.v1);
v0_sol(x,d,Lcz) = simplify(v.v0);
v2_sol(x,d,Lcz) = simplify(v.v2);
phi0_sol(x,d,Lcz) = diff(v0_sol,x);
phi1_sol(x,d,Lcz) = diff(v1_sol,x);
phi2_sol(x,d,Lcz) = diff(v2_sol,x);
M0_sol(x,d,Lcz) = E*I*diff(v0_sol,x,2);
M1_sol(x,d,Lcz) = E*I*diff(v1_sol,x,2);
M2_sol(x,d,Lcz) = E*I*diff(v2_sol,x,2);
T0_sol(x,d,Lcz) = E*I*diff(v0_sol,x,3);
T1_sol(x,d,Lcz) = E*I*diff(v1_sol,x,3);
T2_sol(x,d,Lcz) = E*I*diff(v2_sol,x,3);
%---------
d_lim = solve(v0_sol(0,d,0) == d_0/2,d);
% d_max = solve(v0_sol(0,d,0) == d_0/2,d);
% Lcz_max = solve(v2_sol(-Lcz,d_max,Lcz) - d_c/2 == 0, x,[0 50]);
[d_max, Lcz_max] = solve([v1_sol(0,d,Lcz) - d_0/2 == 0, v2_sol(-Lcz,d_max,Lcz) - d_c/2 == 0],[d,Lcz]);

[Torsten]

Do you succeed with

v = dsolve(odesys);

?

I cannot run your code since variables are not defined.

[EDOARDO GELMI]

This are the costants, and yes i tried without the bc and it works but of course i need to fine the integration costants

E = 70000;           %Young's Modulus [MPa]
h = 6;               %Section Height [mm]
b = 25;              %Section Width [mm]
A = b*h;             %Section Area [mm^2]
I = (b*h^3)/12 ;     %Inertial Momentum [mm^4]
L = 200;             %Specimen Lenght [mm]
ni = 0.33;           %Poisson's Ratio [-]
G = E/(2*(1 + ni));  %Shear Modulus [MPa]
mu = 0.5;            %Shear Factor [-]
A_t = mu*A;          %Shear Section Area [mm^2]
sigma_max = 10;      %Softening Tension [MPa]
d_0 = 0.025;         %Softening Opening [mm]
d_c = 0.1;           %Critical Opening [mm]
a0 = 50;             %Initial Crack Lenght [mm]
%---------
lambda = ((2*b*sigma_max)/(E*I*d_0))^0.25;                      %Constant depending on the bending stiffness of DCB arms and the softening part
k = ((2*b*sigma_max)/(E*I*(d_c-d_0)))^0.25;                     %Costant depending on the linear-elastic stiffness of the interface
w = ((E*I)/(2*G*A_t))*(lambda)^0.5 == 0;                        %Constant depending on bending and shear stiffness of DCB arms, and lambda
ps = ((E*I)/(2*G*A_t))*(k)^0.5 == 0;                            %Constant depending on bending and shear stiffness of DCB arms, and k
%---------
F_lim = (E*I*d_0*lambda^3)/(2*(a0*lambda + sqrt(2*(1+w))));     %Force at which the damage start propagating in the interface
%---------
f1 = linspace(0,F_lim,51);                                      %Generic interval of force in order to simulate the prescribed displacement test in the LE phase 
d1 = ((2*f1*a0)/(3*E*I))*(1+((3*sqrt(2*(1+w)))/(a0*lambda^3))*(a0*lambda*sqrt(2*(1+w))+1+(a0*lambda)^2)); %Crack mouth opening

[Torsten]

In the below code, you have to solve a 12x12 linear system of equations to get the integration constants as functions of d and Lcz. This can take a while (if it is successful at all).

It would be much easier if you could prescribe numerical values for d and Lcz.

E = 70000;           %Young's Modulus [MPa]
h = 6;               %Section Height [mm]
b = 25;              %Section Width [mm]
A = b*h;             %Section Area [mm^2]
I = (b*h^3)/12 ;     %Inertial Momentum [mm^4]
L = 200;             %Specimen Lenght [mm]
ni = 0.33;           %Poisson's Ratio [-]
G = E/(2*(1 + ni));  %Shear Modulus [MPa]
mu = 0.5;            %Shear Factor [-]
A_t = mu*A;          %Shear Section Area [mm^2]
sigma_max = 10;      %Softening Tension [MPa]
d_0 = 0.025;         %Softening Opening [mm]
d_c = 0.1;           %Critical Opening [mm]
a0 = 50;             %Initial Crack Lenght [mm]
%---------
lambda = ((2*b*sigma_max)/(E*I*d_0))^0.25;                      %Constant depending on the bending stiffness of DCB arms and the softening part
k = ((2*b*sigma_max)/(E*I*(d_c-d_0)))^0.25;                     %Costant depending on the linear-elastic stiffness of the interface
w = ((E*I)/(2*G*A_t))*(lambda)^0.5 == 0;                        %Constant depending on bending and shear stiffness of DCB arms, and lambda
ps = ((E*I)/(2*G*A_t))*(k)^0.5 == 0;                            %Constant depending on bending and shear stiffness of DCB arms, and k
%---------
F_lim = (E*I*d_0*lambda^3)/(2*(a0*lambda + sqrt(2*(1+w))));     %Force at which the damage start propagating in the interface
%---------
f1 = linspace(0,F_lim,51);                                      %Generic interval of force in order to simulate the prescribed displacement test in the LE phase 
d1 = ((2*f1*a0)/(3*E*I))*(1+((3*sqrt(2*(1+w)))/(a0*lambda^3))*(a0*lambda*sqrt(2*(1+w))+1+(a0*lambda)^2)); %Crack mouth opening
%% ANALYTICAL MODEL FOR A DCB SPECIMEN UNDER THE CONDITION OF PRESCRIBED DISPLACEMENTS
%% Linear Elastic Phase
%---------
syms x d v0(x) v1(x) v2(x) Lcz 
%---------
ode_0 = diff(v0,x,4) == 0
ode_1 = diff(v1,x,4) - 2*w*(lambda^2)*diff(v1,x,2) + (lambda^4)*v1 == 0
ode_2 = diff(v2,x,4) + 2*ps*(k^2)*diff(v2,x,2) - k^4*(v2 - d_c/2) == 0
sol = dsolve([ode_0,ode_1,ode_2]);
v0(x) = sol.v0;
v1(x) = sol.v1;
v2(x) = sol.v2;
phi0 = -diff(v0,x);
M0 = E*I*diff(v0,x,2);
T0 = E*I*diff(v0,x,3);
phi1 = -diff(v1,x);
M1 = E*I*diff(v1,x,2);
T1 = E*I*diff(v1,x,3);
phi2 = -diff(v2,x);
M2 = E*I*diff(v2,x,2);
T2 = E*I*diff(v2,x,3);

xL = -a0 - Lcz;
xI = -Lcz;
xR = L - a0 - Lcz;
c1 = v0(xL) == d/2;
c2 = M0(xL) == 0;
c3 = v0(xI) == v2(xI);
c4 = phi0(xI) == phi2(xI);
c5 = M0(xI) == M2(xI);
c6 = T0(xI) == T2(xI);
c7 = v1(0) == v2(0);
c8 = phi1(0) == phi2(0);
c9 = M1(0) == M2(0);
c10 = T1(0) == T2(0);
c11 = v1(xR) == 0;
c12 = phi1(xR) == 0;
c = [c1,c2,c3,c4,c5,c6,c7,c8,c9,c10,c11,c12];
vars = symvar(c);
[A,b] = equationsToMatrix(c,'Vars',vars(1:12));
constants = A\b
%constants = solve(c,v(1:12));   % very time-consuming
v1(x) = subs(v1(x),vars(1:12),constants)
v2(x) = subs(v2(x),vars(1:12),constants)
v3(x) = subs(v3(x),vars(1:12),constants)

[EDOARDO GELMI]

Yeah i know it should work faster with a numerical solution but unfortunatly i cannot use it. I can try your script, right now mine is working but it's very time consuming (it's been an hour untill now)

[EDOARDO GELMI]

Thank you very much, i'll let you know if it works