ploting second derivative from second order differential equationn

2021 3 18
1 个回答
回答已采纳
更新时间：2021 3 20
3 次查看（30 天）

0 个投票

Hi, Beloow is my code and i am solving second order differential equation numericaly and i am getting plot of first and zeroth derivation from ode45. is there any way to plot the second derivative as well by just calling ode45 ?
function fval = simulation( t,y )
x=y(1);
v=y(2);
c=6850;
m=450;
k=0;
f=2650;
e=535;
w=4410;
fval(1,1)=v;
%fval(2,1)= -(c*v+k*x)/m
fval(2,1)=(-c*v-f-e+w+k*x)/m;
end
%%%%main file%%%%%%
c=6850;
m=450;
k=0;
f=2650;
e=535;
w=4410;
Y0= [0;3.86];
tspan= [0 1];
options= odeset ('Events', @event);
[tsol,ysol] = ode45(@(t,y) simulation(t,y), tspan, Y0,options);
figure
plot(tsol,ysol(:,1), 'r-');
hold on;
figure
plot(tsol,ysol(:,2), 'g-');
hold on;
figure

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

 采纳的回答

Jan
Jan 2021-3-18

1 个投票

The easiest solution is to let the function to be integrated accept vectors as input and call it with the time steps and positions replied by ODE45:
[tsol, ysol] = ode45(@(t,y) simulation(t,y), tspan, Y0,options);
dy = simulation(tsol, ysol.').';
d2y = dy(:, 2);
function fval = simulation(t, y)
x = y(1, :);
v = y(2, :);
c = 6850;
m = 450;
k = 0;
f = 2650;
e = 535;
w = 4410;
fval(1, :) = v;
fval(2, :) = (-c * v - f - e + w + k * x) / m;
end

12 个评论
显示 10更早的评论 隐藏 10更早的评论

Jan there is some syntax error in your code line 2
Jan
Jan 2021-3-18
编辑：Jan 2021-3-18
The 2nd line of my code is empty. So please post the error message, if you get one. Letting me guess which error appears anywhere in the code, is not efficient.
By the way: My code contains only the snippet you have to change to get the 2nd derivative. You have to insert this in your code.
Muhammad Sarmad Zahid
编辑：Jan 2021-3-19
i have one more question i have this code. I am following you guidance to plot g against time. g is a contant value. How can i plot g ?
Y0= [-0.75 ;0];
tspan= [0 2];
options= odeset ('Events', @freefall);
[tsol,ysol] = ode45(@(t,y) pra(t,y), tspan, Y0,options);
figure
plot(tsol,ysol(:,1), 'r-');
hold on;
figure
plot(tsol,ysol(:,2), 'g-');
hold on;
figure
%%%%function%%%
function fval = pra( t,y )
g=9.8;
y(1);
y(2);
% c=6850;
fval = [y(2); g)];
end
Jan
Jan 2021-3-19
编辑：Jan 2021-3-19
There is a not matching parenthesis in your code in:
fval = [y(2); g)];
You can plot the 2nd component of the output of pra() in exactly the same way, that I've showed you for simulation(): Convert pra() such that it accept vectors as input and use the output of ODE45 to get the wanted values. Afterwards use plot() as for the other components.
Muhammad Sarmad Zahid
编辑：Jan 2021-3-19
i am tring to do that but its not returning me a contant value of g instead it is returning me big values starting from 2.7 to 55 m/s^2. can you please guideme by writing a code? it would be really appreaciat if you check that code as well.
you can chekc the values of f as well
[t,y] = ode45(@(t,y) pra(t,y), [0 1], [-0.75 ; 0 ]);
f=simulation( t,y.').'
%%%%function%%%
function fval = pra( t,y )
g=9.8;
y(1);
y(2);
% c=6850;
fval = [y(2); g.*((t).^0)];
end
Jan
Jan 2021-3-19
编辑：Jan 2021-3-19
Please use the tools to format the code. This makes it easier to read you code. I've done this for you, but it is more efficient, if you do this by your own.
I've showed you how to convert the function simulation() such, that it accepts matrices as input. Then you can use the output of ODE45 to obtain the 2nd derivative also.
Now you have to modify your function pra() accordingly.
These lines are meaningless:
y(1);
y(2);
Instead of:
fval = [y(2); g.*((t).^0)];
you need:
fval = [y(2, :); repmat(g, size(t))];
Using g.*((t).^0) is smart, but the power operation is very expensive. repmat() is more efficient, or alternatively:
fval = [y(2, :); g * ones(size(t))];
Afterwards, do not call "f=simulation( t,y.').'", but:
df = pra(t, y.').';
plot(t, df(:, 2))
Error using vertcat
Dimensions of arrays being concatenated are not consistent.
Error in pra (line 10)
fval = [y(2, :); g * ones(size(t))];
Error in practicemailfile (line 20)
df = pra(tsol, ysol.').';
Jan
Jan 2021-3-19
Then call it like this:
df = pra(tsol.', ysol.').';
% ^^
Jan
Jan 2021-3-20
But this is the same question again?

请先登录，再进行评论。

更多回答（0 个）

类别

在 帮助中心 和 File Exchange 中查找有关 Solver Outputs and Iterative Display 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by