Second derivative of a function handle

Hi, I have a function that I want to differentiate twice and plot it against time for different values of a. I looked through some forum posts in which they recommend to use syms but I haven't learned about syms yet. I tried the following code but the system will return an error of "This statement is incomplete." Can someone have a look?

syms t
D = 5; % Damping factor
k = 37; % Spring constant
m = 0.5; % Quarter mass of the car
H = 0.1; % Height
L = 0.2; % Length
t = (0:1:50);
for a = [1, 5, 25, 50, 100] 
    b = (-(k - m*a^2)*k*H)/((k - m*a^2)^2 + (a*D)^2);
    c = (-a*D*k*H)/((k-m*a^2)^2 + (a*D)^2);
    zp = @(t)b*cos(a.*t) + c*sin(a.*t) + H;
    zp1 = eval(['@(t)' char(diff(zp(t)))]);
    zp2 = eval(['@(t)' char(diff(zp1(t)))]);
    plot(t,zp2(t));
    hold on 
end

[Stanley Ka]

Thanks a lot, I'll have a go with what you said.

Answer 1

[Stephen23]

Tip for beginners: whenever you use EVAL you are doing something wrong.

With symbolic variables you do not need to define anonymous functions like that.

This should get you started, you probably need to fix the code in other ways (note that you will not get anything meaningful out of that plot when the data varies so much in amplitude). I recommend storing the data in a matrix and plotting after the loop (makes it easier to add a legend).

syms t
D = 5; % Damping factor
k = 37; % Spring constant
m = 0.5; % Quarter mass of the car
H = 0.1; % Height
L = 0.2; % Length
T = 0:0.1:50;
for a = [10,5,1] 
    b = (-(k - m*a^2)*k*H)/((k - m*a^2)^2 + (a*D)^2);
    c = (-a*D*k*H)/((k-m*a^2)^2 + (a*D)^2);
    zp0 = b*cos(a.*t) + c*sin(a.*t) + H;
    zp1 = diff(zp0);
    zp2 = diff(zp1);
    Z = double(subs(zp2,t,T));
    plot(T,Z)
    hold on
end

Answer 2

[Sulaymon Eshkabilov]

Hi,

Here is a corrected and more efficient solution to your exercise.

clc
clearvars
syms zp(t) a b c
D = 5; % Damping factor
k = 37; % Spring constant
m = 0.5; % Quarter mass of the car
H = 0.1; % Height
L = 0.2; % Length
zp(t) = b*cos(a.*t) + c*sin(a.*t) + H;
Dzp=diff(zp(t),t);
D2zp=diff(Dzp,t);
t = (0:.001:1); % t = [0:1:50] is not a good resolution and thus, t = [0, 1] is chosen
%%
a =double([1, 5, 25, 50, 100]);
b = (-(k - m*a.^2)*k*H)./((k - m*a.^2).^2 + (a*D).^2);
c = (-a*D*k*H)./((k-m*a.^2).^2 + (a*D).^2);
DC=(- a(:).^2.*b(:).*cos(a(:).*t) - (a(:).^2).*c(:).*sin(a(:).*t)); % This is the 2nd diff of zp(t)
plot(t, DC, '-')

[Torsten]

And why do you calculate the second derivative using the diff command if finally you use a computation by hand ?

[Sulaymon Eshkabilov]

That is just copy and paste of the computed formulation from diff() above.