# why do I get an error "First input argument must be a function handle"

1 view (last 30 days)
Cedrick De Block on 23 Nov 2021
Answered: Star Strider on 23 Nov 2021
Hello,
I want to calculate the integral of the function F_z, like this:
nUn = 2000;
nU2n = 1500;
nU3n = 200;
nU4n = 20;
nVn = nUn;
nV2n = nU2n;
nV3n = nU3n;
nV4n = nU4n;
dt = 0.001 ; % Timestep
t = 0:dt:1; % Time interval of the vertical excitation
Om = 405 * 2 * pi / 60; % Frequency of exitation [rpm]
phi = Om * t;
n = 3; % Number of blades
F_z = -(nUn*cos(n*phi) + nU2n*cos(2*n*phi) + nU3n*cos(3*n*phi) + ...
+ nU4n*cos(4*n*phi) + nVn*sin(n*phi) + nV2n*sin(2*n*phi) + ...
+ nV3n*sin(3*n*phi) + nV4n*sin(4*n*phi));
% Plot of the evolution fo the vertical excitation
plot(t, F_z)
% Computing and plot of PSD
for f = 0:0.1:100
F_z_sqMean = f * integral((F_z).^2, 0, 1/f);
end
##### 0 CommentsShowHide -1 older comments

Sign in to comment.

### Answers (1)

Star Strider on 23 Nov 2021
I have no idea what the code does.
See if this paproximates the desired result —
nUn = 2000;
nU2n = 1500;
nU3n = 200;
nU4n = 20;
nVn = nUn;
nV2n = nU2n;
nV3n = nU3n;
nV4n = nU4n;
dt = 0.001 ; % Timestep
t = 0:dt:1; % Time interval of the vertical excitation
Om = 405 * 2 * pi / 60; % Frequency of exitation [rpm]
phi = @(t) Om * t; % <— Is This Part Of The Desired Outcome?
n = 3; % Number of blades
F_z = @(t) -(nUn*cos(n*phi(t)) + nU2n*cos(2*n*phi(t)) + nU3n*cos(3*n*phi(t)) + ...
+ nU4n*cos(4*n*phi(t)) + nVn*sin(n*phi(t)) + nV2n*sin(2*n*phi(t)) + ...
+ nV3n*sin(3*n*phi(t)) + nV4n*sin(4*n*phi(t)));
% Plot of the evolution fo the vertical excitation
figure
plot(t, F_z(t))
grid
xlabel('t')
ylabel('F_z(t)')
% Computing and plot of PSD
f = 0:0.1:100;
for k = 1:numel(f)
F_z_sqMean(k,:) = f * integral(@(t)F_z(t).^2, 0, 1/f(k));
end
Warning: Reached the limit on the maximum number of intervals in use. Approximate bound on error is 1.1e+33. The integral may not exist, or it may be difficult to approximate numerically to the requested accuracy.
figure
surfc(F_z_sqMean, 'EdgeColor','none')
grid on
set(gca, 'XScale','log', 'YScale','log', 'ZScale','log')
view(115,30)
xlabel('X')
ylabel('Y')
Experiment!
.
##### 0 CommentsShowHide -1 older comments

Sign in to comment.

### Community Treasure Hunt

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

Start Hunting!