Solve differential equation to get variable

% I am trying to proceed to get the solution for a here by equating differentiation of u_a withrespect to a equal to Zero, but got error like!

%Warning: Solutions are only valid under certain conditions. To include parameters and %conditions in the solution, specify the 'ReturnConditions' value as

%'true'.

%> In sym/solve>warnIfParams (line 478)

%In sym/solve (line 357)

%In HW2_Pro2_VirtualForce (line 55)

%>>

clc;

close;

clear;

syms d_F ;

syms d_F_c;

syms dM_B;

%syms Theta_B;

syms x real;

syms u_mid;

syms u_c;

syms u_max;

syms a;

%syms w_x;

syms M_z(x);

syms dM_z(x);

syms dP;

syms u_a;

%syms n;

%Q = 2.5; %kip/ft

L = 30.; % ft

%a = 18; %assume maximum deflection at 15', 20'

E = 29000.*12^2; %ksi to kips/ft2

I = 1890./(12)^4; %ft^4

E_I = E.*I; %kip*ft^2;

%% deflection at mid of AB

w_x = @(x) 1*x/(L/3)*heaviside(x-L/3)*(1-heaviside(x-2*L/3));

w_x_1 = w_x(x);

eq1 = diff(M_z(x),2) == -w_x_1;

sol1(x) = dsolve(eq1, M_z(0)==-150, M_z(L) == 250, 'IgnoreAnalyticConstraints',true);

set(gca, 'XAxisLocation', 'origin', 'YAxisLocation', 'origin')

fplot(x, sol1(x), [0,L]);

figure();

%moment_x = M_z(x)/sol1;

fplot(x, w_x(x), [0,L]);

origin = [0 0];

figure();

%%

%apply virtual load at distance a from support A

%a = L/2;

%dP = -1.;

dM_z= @(x,a) (dP*x*(L-a)/L)-dP.*(x-a).*heaviside(x-a);

%fplot(x, dM_z(x,a), [0,L])

dW = @(x,a) M_z(x)*dM_z(x,a);

dW_int = @(a) 1/E_I.*(int(dW, x, 0 ,L));

dW_ext = dP*u_a;

eqn = dW_int == dW_ext;

sol = solve(eqn, u_a);

u_a_diff = diff(u_a, a);

eqn2 = u_a_diff == 0;

sol_a = solve(eqn2, a);

a_solve = double(sol_a);

Answer 1

[Paul]

Hi Milan,

I modified as follows. Not sure if this is what you're looking for:

syms d_F ;
syms d_F_c;
syms dM_B;
%syms Theta_B;
syms x real;
syms u_mid;
syms u_c;
syms u_max;
syms a;
%syms w_x;
syms M_z(x);
syms dM_z(x);
syms dP;
syms u_a;
%syms n;
%Q = 2.5; %kip/ft
L  = 30.; % ft 
%a = 18;  %assume maximum deflection at 15', 20' 
E = 29000.*12^2; %ksi to kips/ft2
I = 1890./(12)^4; %ft^4
E_I = E.*I; %kip*ft^2;
%% deflection at mid of AB
w_x = @(x) 1*x/(L/3)*heaviside(x-L/3)*(1-heaviside(x-2*L/3));
w_x_1 = w_x(x);
eq1 = diff(M_z(x),2) == -w_x_1;
sol1(x) = dsolve(eq1, M_z(0)==-150, M_z(L) == 250, 'IgnoreAnalyticConstraints',true)

sol1(x) = 

sol1(x) is the solution for M_z(x), so make that assignment because M_z(x) is used below.

M_z(x) = sol1(x);

% set(gca, 'XAxisLocation', 'origin', 'YAxisLocation', 'origin')
% fplot(x, sol1(x), [0,L]);
% figure();
%moment_x = M_z(x)/sol1;
% fplot(x, w_x(x), [0,L]);
% origin = [0 0];
% figure();
%% 
%apply virtual load at distance a from support A
%a  = L/2;
%dP = -1.;

Use symfun objects in instead of anonymous function

%dM_z= @(x,a) (dP*x*(L-a)/L)-dP.*(x-a).*heaviside(x-a);
dM_z(x,a) = (dP*x*(L-a)/L)-dP.*(x-a).*heaviside(x-a)

dM_z(x, a) = 

%fplot(x, dM_z(x,a), [0,L])

%dW = @(x,a) M_z(x)*dM_z(x,a);
dW(x,a) = M_z(x)*dM_z(x,a);
 
dW_int(a) =  1/E_I.*(int(dW(x,a), x, 0 ,L));
dW_ext = dP*u_a; 
eqn = dW_int == dW_ext;

Solve eqn for u_a

sol = solve(eqn, u_a, 'ReturnConditions',true)

sol = struct with fields:
           u_a: [6×1 sym]
    parameters: [1×0 sym]
    conditions: [6×1 sym]

We see six possible solutions depending on the value of a, all with dP ~= 0.

sol.conditions

ans = 

Assume that a is bracketed by 10 < a < 20 because there was a comment upstream that a = 18. Assume that dP ~= 0

assumeAlso(10 < a < 20);
assumeAlso(dP ~= 0);
sol = solve(eqn, u_a, 'ReturnConditions',true);

Under those assumptions, we have one solution for u_a

sol.u_a

ans = 

Differentiate u_a wrt a

u_a_diff = diff(sol.u_a, a)

u_a_diff = 

And solve for a

eqn2 = u_a_diff == 0;
sol_a = solve(eqn2, a);
a_solve = double(sol_a)

a_solve = 18.1398

[Milan]

Thanks alot Paul this worked out. However, I tried to get value of u_a at a_solve, but got an error Error using sym/subsindex

Invalid indexing or function definition. Indexing must follow MATLAB indexing.

Function arguments must be symbolic variables, and function body must be sym

expression.

a_solve = double(sol_a)

M_max = double(M_z(a_solve))

u_max = u_a_1(a_solve)

[Paul]

Hard to know what's going on w/o seeing the actual code that generated the error message. See below for what I think the goal is. Not that u_a_1 is not defined in the original code, so I had to guess what you're trying to do.

syms d_F ;
syms d_F_c;
syms dM_B;
%syms Theta_B;
syms x real;
syms u_mid;
syms u_c;
syms u_max;
syms a;
%syms w_x;
syms M_z(x);
syms dM_z(x);
syms dP;
syms u_a;
%syms n;
%Q = 2.5; %kip/ft
L  = 30.; % ft 
%a = 18;  %assume maximum deflection at 15', 20' 
E = 29000.*12^2; %ksi to kips/ft2
I = 1890./(12)^4; %ft^4
E_I = E.*I; %kip*ft^2;
%% deflection at mid of AB
w_x = @(x) 1*x/(L/3)*heaviside(x-L/3)*(1-heaviside(x-2*L/3));
w_x_1 = w_x(x);
eq1 = diff(M_z(x),2) == -w_x_1;
sol1(x) = dsolve(eq1, M_z(0)==-150, M_z(L) == 250, 'IgnoreAnalyticConstraints',true);

M_z(x) = sol1(x);

dM_z(x,a) = (dP*x*(L-a)/L)-dP.*(x-a).*heaviside(x-a);

dW(x,a) = M_z(x)*dM_z(x,a);
 
dW_int(a) =  1/E_I.*(int(dW(x,a), x, 0 ,L));
dW_ext = dP*u_a; 
eqn = dW_int == dW_ext;

assumeAlso(10 < a < 20);
assumeAlso(dP ~= 0);
sol = solve(eqn, u_a, 'ReturnConditions',true);

Define u_a

u_a(a) = sol.u_a;

u_a_diff = diff(u_a, a);

eqn2 = u_a_diff == 0;
a_solve = solve(eqn2, a);
a_solve = double(a_solve)

a_solve = 18.1398

M_max = double(M_z(a_solve))

M_max = 180.7568

%u_max = u_a_1(a_solve) % u_a_1 is not defined!
u_max = double(u_a(a_solve))

u_max = 0.0380

[Milan]

This worked Paul Thanks alot for your help, I appreciate it!