getting output in the form of 'vpaintegral' when applying dsolve command

3 visualizaciones (últimos 30 días)
aakash dewangan
aakash dewangan el 29 de Mzo. de 2022
Respondida: Walter Roberson el 11 de Abr. de 2022
Hello All,
I have written a code to solve a problem. In my code, in the last section, when I am applying "dsolve" command, code gives me an error. Could anybody please help me to solve this error.
Waiting for responses.
Thank you
%% AAKASH DEWANGAN 9/12/2021
clc; clear all; close all
syms p1(t) p2(t) p3(t) p4(t) p5(t) p6(t) rho L m v T k G y_mass(t)
%% parameters
rho = 1.3; T = 45000; L = 60; k = 15000; m = 7; v = 350*1000/3600; G = 0.1 % HIGH speed parameters
Dp1 = diff(p1); D2p1 = diff(p1,2); Dp2 = diff(p2); D2p2 = diff(p2,2); Dp3 = diff(p3); D2p3 = diff(p3,2);
%%
% first matrix terms
AA = rho*L/2 + m*(sin(pi*v*t/L))^2;
BB = m*sin(2*pi*v*t/L)*sin(pi*v*t/L);
CC = m*sin(3*pi*v*t/L)*sin(pi*v*t/L);
DD = rho*L/2 + m*(sin(2*pi*v*t/L))^2;
EE = m*sin(2*pi*v*t/L)*sin(3*pi*v*t/L);
FF = rho*L/2 + m*(sin(3*pi*v*t/L))^2;
% second matrix terms
GG = T*(pi/L)^2*(L/2) + k*(sin(pi*v*t/L))^2;
HH = k*sin(2*pi*v*t/L)*sin(pi*v*t/L);
II = k*sin(pi*v*t/L)*sin(3*pi*v*t/L);
JJ = T*(2*pi/L)^2*(L/2) + k*(sin(2*pi*v*t/L))^2;
KK = k*sin(2*pi*v*t/L)*sin(3*pi*v*t/L);
LL = T*(3*pi/L)^2*(L/2) + k*(sin(3*pi*v*t/L))^2;
% RHS
MM = k*G*sin(pi*v*t/L);
NN = k*G*sin(2*pi*v*t/L);
OO = k*G*sin(3*pi*v*t/L);
%%
tim = zeros(1,2)
poly_order = 10;
F_val = k*G; jloss = 1;
y_massVal = G
p3_expression = 0; p1_expression = 0; p2_expression = 0; p3dot_expression = 0; p1dot_expression = 0; p2dot_expression = 0;
axis tight
vid = VideoWriter('ForMATLAB.avi');
open(vid);
path = pwd ;
%%
for contacts = 1:100
contacts
% Equation (coupled system of ODE to solve for p)
Eq1 = AA*diff(p1,t,2) + BB*diff(p2,t,2) + CC*diff(p3,t,2) + GG*p1 + HH*p2 + II*p3 == MM; % Equation 1
Eq2 = BB*diff(p1,t,2) + DD*diff(p2,t,2) + EE*diff(p3,t,2) + HH*p1 + JJ*p2 + KK*p3 == NN; % Equation 2
Eq3 = CC*diff(p1,t,2) + EE*diff(p2,t,2) + FF*diff(p3,t,2) + II*p1 + KK*p2 + LL*p3 == OO; % Equation 3
%%
[V,S] = odeToVectorField(Eq1, Eq2, Eq3); % converts ODE in state space form
ftotal = matlabFunction(V, 'Vars',{'t','Y'}); % Using readymade MATLAB function to solve using ODE 45
% ^-^ - single quotes1 + l*p2 + m*p3== p
tstart = tim(jloss)
interval = [tstart L/v]; % Time Interval to run the program
p3_IC = subs(p3_expression,t,tim(jloss))
p3dot_IC = subs(p3dot_expression,t,tim(jloss))
p2_IC = subs(p2_expression,t,tim(jloss))
p2dot_IC = subs(p2dot_expression,t,tim(jloss))
p1_IC = subs(p1_expression,t,tim(jloss))
p1dot_IC = subs(p1dot_expression,t,tim(jloss))
IC = double([p3_IC p3dot_IC p1_IC p1dot_IC p2_IC p2dot_IC ])
%% ==========================================================================
[tim pSol] = ode45(@(t,Y)ftotal(t,Y),interval,IC); % Using ODE 45 to solve stste space form of ODE
p3Values = (pSol(:,1)); % number 1 denotes first solution likewise you can mention 2 ,3 & 4 for the next three solutions
p3dotValues = (pSol(:,2));
p1Values = (pSol(:,3)); % number 1 denotes first solution likewise you can mention 2 ,3 & 4 for the next three solutions
p1dotValues = (pSol(:,4));
p2Values = (pSol(:,5)); % number 1 denotes first solution likewise you can mention 2 ,3 & 4 for the next three solutions
p2dotValues = (pSol(:,6));
%% Curve fitting
p_1 = polyfit(tim,p1Values,poly_order) % curve fitting of data points using polynomial (third argument shows degree of polynomial)
for i = 1:length(p_1)
ele_p_1(i) = p_1(i)*t^(length(p_1)-i);
end
p1_expression = vpa(sum(ele_p_1));
%%
p_1dot = polyfit(tim,p1dotValues,poly_order) % curve fitting of data points using polynomial (third argument shows degree of polynomial)
for idot = 1:length(p_1dot)
ele_p_1dot(idot) = p_1dot(idot)*t^(length(p_1dot)-idot);
end
p1dot_expression = vpa(sum(ele_p_1dot));
p_2 = polyfit(tim,p2Values,poly_order) % curve fitting of data points using polynomial (third argument shows degree of polynomial)
for j = 1:length(p_2)
ele_p_2(j) = p_2(j)*t^(length(p_2)-j);
end
p2_expression = vpa(sum(ele_p_2));
%%
p_2dot = polyfit(tim,p2dotValues,poly_order); % curve fitting of data points using polynomial (third argument shows degree of polynomial)
for jdot = 1:length(p_2dot)
ele_p_2dot(jdot) = p_2dot(jdot)*t^(length(p_2dot)-jdot);
end
p2dot_expression = vpa(sum(ele_p_2dot));
p_3 = polyfit(tim,p3Values,poly_order) % curve fitting of data points using polynomial (third argument shows degree of polynomial)
for ii = 1:length(p_3)
ele_p_3(ii) = p_3(ii)*t^(length(p_3)-ii);
end
p3_expression = vpa(sum(ele_p_3));
p_3dot = polyfit(tim,p3dotValues,poly_order) % curve fitting of data points using polynomial (third argument shows degree of polynomial)
for iidot = 1:length(p_3dot)
ele_p_3dot(iidot) = p_3dot(iidot)*t^(length(p_3dot)-iidot);
end
p3dot_expression = vpa(sum(ele_p_3dot));
%% Displacement u
syms x
ter1(t) = sin(pi*x/L)*p1_expression;
ter2(t) = sin(2*pi*x/L)*p2_expression;
ter3(t) = sin(3*pi*x/L)*p3_expression;
u = ter1 + ter2 + ter3
%% Force
u_vt = subs(u,x,v*t);
dd_u_vt = diff(u_vt,t,2);
F = k*G-k*u_vt-m*dd_u_vt
break
end
%% Error part (this part is giving me error)
syms y(t)
Dy = diff(y,t)
equation = m*diff(y,t,2) + k*y == F
condtn1 = y(tstart) == G; condtn2 = Dy(tstart) == 0;
condition = [condtn1; condtn2]
sol_y_mass = dsolve(equation,condition)
you = G - sol_y_mass
figure(101)
ezplot(you,[tstart L/v])
hold on
xlim([0,tim(L/v)])
beep
  13 comentarios
Mostrar 11 comentarios más antiguos
Torsten
Torsten el 8 de Abr. de 2022
Does
sol_y_mass = matlabFunction(sol_y_mass)
work ?
If yes, are there other variables except t that appear in the function handle ?
aakash dewangan
aakash dewangan el 11 de Abr. de 2022
No, It did Not work.
There are no other variables except t in the function.
The best way to help is to run the code in your computer, if possible. :)
Thanks

Iniciar sesión para comentar.

Iniciar sesión para responder a esta pregunta.

Respuestas (1)

Walter Roberson
Walter Roberson el 11 de Abr. de 2022
Change the plotting to something like this:
tvec = linspace(tstart, L/v, 200);
Y = double(subs(you, t, tvec));
plot(tvec, Y)
hold on
xlim([0,(L/v)])
You will definitely not be able to get ezplot() to work.
You can try fplot() instead of ezplot(), but it would take rather a long time.

Categorías

Más información sobre Equation Solving en Help Center y File Exchange.

Etiquetas

Productos


Versión

R2021a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by