Solve differential equation to get variable
1 view (last 30 days)
Show older comments
% 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);
0 Comments
Accepted Answer
Paul
on 30 Oct 2022
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) 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)
%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)
We see six possible solutions depending on the value of a, all with dP ~= 0.
sol.conditions
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
Differentiate u_a wrt a
u_a_diff = diff(sol.u_a, a)
And solve for a
eqn2 = u_a_diff == 0;
sol_a = solve(eqn2, a);
a_solve = double(sol_a)
3 Comments
Paul
on 30 Oct 2022
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)
M_max = double(M_z(a_solve))
%u_max = u_a_1(a_solve) % u_a_1 is not defined!
u_max = double(u_a(a_solve))
More Answers (0)
See Also
Categories
Find more on Equation Solving in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!