Solve differential equation to get variable
1 Ansicht (letzte 30 Tage)
Ältere Kommentare anzeigen
% 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 Kommentare
Akzeptierte Antwort
Paul
am 30 Okt. 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 Kommentare
Paul
am 30 Okt. 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))
Weitere Antworten (0)
Siehe auch
Kategorien
Mehr zu Equation Solving finden Sie in Help Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!