MATLAB: How to create a trapezoidal rule routine to evaluate the following Bessel integral without Matlab functions. ME 1101

besselMATLABnumerical integrationtrapezoidal

Best Answer

clear all
close all
clc
% Integrates & plots Bessel functions of the first kind
a = 0; % x intial
b = 20;  % x final
c = 0;  % theta intial
d = pi;  % theta final
n = 1;  % order of the Bessel integral
int = 100;  % intervals for x and theta
theta = linspace(c,d,101);
x = linspace(a,b,101);
for i = 1:length(x)
    J_x = (1/pi)*cos(n*theta-x(i)*sin(theta));
    J(i) = (d-c)*(sum(J_x)+sum(J_x(2:length(J_x)-1)))/(2*int);
   
end
plot(x,J)
grid on
grid minor
title('Bessel Function of the First Kind')
xlabel('x')
ylabel('Jn(x)')
% Check with the Bessel function
y = besselj(1,x);
figure
plot(x,y,'b')
grid on
grid minor
title('Bessel Function of the First Kind')
xlabel('x')
ylabel('J_{1}(x)')

Related Solutions

MATLAB: Problem in simulation ODE45 Index exceeds array bounds.

 I added more elements to the code and also new conditions but this one worked for me.
    % Function
function [xdot y] = HW5P3(t,x)% x is theta
%Function equations
global a b B J K vo R Lo thetaout
theta=x(1);
h=x(2);
lambda=x(3);
v=(t<0.5).*(0)+((t>=0.5)&(t<4.5)).*(vo*sin(2*pi*(1/0.1)*t))+(t>=4.5).*(0); %Voltage source starts at t=0.5sec
L=(theta<-thetaout).*(Lo)+((theta>=-thetaout)&(theta<thetaout)).*(a+b.*abs(theta))+(theta>=thetaout).*(Lo); %Inductance change with theta
Ldot=b*sign(theta);%Derivative of L
lambdadot=v-R*(lambda/L);
hdot=((Ldot/(2*(L)^2))*lambda^2)-((B/J)*h)-K*theta;
thetadot=h/J;
xdot=[thetadot;hdot;lambdadot];
y=[L;v;Ldot];
end
      % code
global a b B J K vo R Lo thetaout
a=1;
b=-2/pi;
B=0.01;%Damping factor for bearings
J=0.001;%Inertia of the shaft
K=1;%Rotational Spring of the shaft
Lo=0.1;%minimum inductance when the iron of the blade is outside the gap
vo=120;%source voltage pk-pk
R=10; %The resistance value was to high and the source wasn't giving enough current. I changed it to 10ohms
thetaout=1.4137;%limit of the blade to work by the effect of L(theta) out of this angle L=Lo
lambdai=0;%Initial condition
thetai=1.4137;%Initial COndition

hi=0;%Initial COndition
[t,x]=ode45('HW5P3',[0 5],[thetai,hi,lambdai]);  %calling of function
L=length(t);
N = length(t);
for i=1:N %getting values of output from the function for every time point
[xd,y] = HW5P3(t(i),x(i,:));
Ln(i) = y(1);%L(theta)
vs(i) = y(2);%Voltage source
Lndot(i) = y(3);%Derivative of L(theta)
end
theta=x(:,1);
h=x(:,2);
lambda=x(:,3);
thetadot=h/J;
is=lambda/L;
Trotspring=K*theta;
Tbearing=B*thetadot;
figure(1);
subplot(3,1,1)
plot(t,theta), xlabel('time[sec]'),ylabel('Angular displacement [rad]')
subplot(3,1,2)
plot(t,thetadot), xlabel('time[sec]'),ylabel('Angular velocity [rad/sec]')
subplot(3,1,3)
plot(t,Ln,t,theta), xlabel('time[sec]'),legend('Ln [H]','Angular displacement [rad]')
figure(2);
subplot(3,1,1)
plot(theta,Ln), xlabel('Angular displacement [rad]'),ylabel('L [H]')
subplot(3,1,2)
plot(t,theta,t,is*1000), xlabel('time[sec]'),legend('Angular displacement [rad]','Current [mA]')
subplot(3,1,3)
plot(t,Trotspring,t,Tbearing), xlabel('time[sec]'),legend('Spring Torque [N]','Bearings Tor

MATLAB: SUBPLOT : UPPER & LOWER PLOT

%true
for i = 1: 6
 subplot(6,2,i)
 y = cos(i*x);
 plot(x,y)
 title(sprintf('y = cos(%dx)', i));
end

Best Answer

Related Solutions

MATLAB: Problem in simulation ODE45 Index exceeds array bounds.

MATLAB: SUBPLOT : UPPER & LOWER PLOT

Related Question