getting output in the form of 'vpaintegral' when applying dsolve command
Show older comments
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
Torsten
on 29 Mar 2022
Does the form of "F" promise success for the integration ?
Could you show us "F" before the call to "dsolve" for "sol_y_mass" ?
aakash dewangan
on 30 Mar 2022
Edited: aakash dewangan
on 30 Mar 2022
Bjorn Gustavsson
on 30 Mar 2022
Can you simply integrate the non-homogenous ODE with a general 10-degree polynomial on the RHS and then identify and tranform/replace the polynomial coefficients to the ones you have at the end of the pre-dsolve step in the script. Since it is a linear ODE with a comparatively simple non-homogenous part this should be simple - on the other hand I think many things "should be simple" that aren't.
Torsten
on 30 Mar 2022
What error message do you get ? In principle, it should be no problem to integrate a 2nd order linear ODE with polynomial right-hand side:
syms m k t y(t)
d2ydt2 = diff(y,t,2);
F = 1 + 2*t^2 + 5*t^3;
sol = dsolve(m*d2ydt2 + k*y == F)
aakash dewangan
on 30 Mar 2022
Torsten
on 30 Mar 2022
I can't run your code because I don't have MATLAB available at the moment.
Could you show the error message and the values for G, m, k and F ?
aakash dewangan
on 3 Apr 2022
Torsten
on 3 Apr 2022
Seems the error comes from the "ezplot" command, not from "dsolve".
aakash dewangan
on 6 Apr 2022
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
youf = matlabFunction(you)
figure(101)
ezplot(youf,[tstart L/v])
hold on
xlim([0,tim(L/v)])
beep
Does this work ?
If not, please show the error message.
aakash dewangan
on 8 Apr 2022
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
on 11 Apr 2022
Answers (1)
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 Mathematics and Optimization in Help Center and File Exchange
on 29 Mar 2022
on 11 Apr 2022
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
