清空筛选条件
清空筛选条件

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

2 次查看（过去 30 天）
aakash dewangan
aakash dewangan 2022-3-29
回答： Walter Roberson 2022-4-11
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 个评论
显示 11更早的评论
Torsten
Torsten 2022-4-8
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 2022-4-11
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

请先登录，再进行评论。

请先登录，再回答此问题。

回答（1 个）

Walter Roberson
Walter Roberson 2022-4-11
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.

类别

在 Help Center 和 File Exchange 中查找有关 Equation Solving 的更多信息

标签

产品


版本

R2021a

Community Treasure Hunt

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

Start Hunting!

Translated by