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

1 view (last 30 days)
aakash dewangan
aakash dewangan on 29 Mar 2022
Answered: Walter Roberson on 11 Apr 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 Comments
Show 11 older comments
Torsten
Torsten on 8 Apr 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 on 11 Apr 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

Sign in to comment.

Sign in to answer this question.

Answers (1)

Walter Roberson
Walter Roberson on 11 Apr 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.

Categories

Find more on Nonlinear Regression in Help Center and File Exchange

Tags

Products


Release

R2021a

Community Treasure Hunt

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

Start Hunting!