Download pendulum.m(at the bottom)(the nonlinear pendulum equation solved via th

Question

Download pendulum.m(at the bottom)(the nonlinear pendulum equation solved via the Euler Method).

1. Copy the original code to create a new script called pendulum_pc.m. Using notes from class, convert this to use the Predictor-Corrector Method (Euler step to estimate derivative at next step, then use the average of this estimate with the Euler value to actually advance). Reduce the stepsize of this method by orders of magnitude (i.e. reduce from 1e-4 to 1e-3, etc) until you see divergence. Is this method more stable than Euler alone?

2. Copy the original code to create a new script called pendulum_mp.m. Convert it now to use the Midpoint Method (MP). In this method, you take a half step (h/2) to the midpoint of the interval between the current step and the next step to estimate the derivative there, then use that derivative alone to advance to the next step. So: f(i+1/2) = f(i) + h/2 * df/dt(i), then f(i+1) = f(i) + h * df/dt(i+1/2). Again here, explore larger stepsizes with respect to stability. Is this method as stable as the Predictor-Corrector Method?

3. Download pend_driver.m and pend.m.(located at the bottom) These of course currently solve the nonlinear pendulum equation using a canned package (ode23). Convert them to solve the coupled linear second order forced system:

Solve this for k1 = 2 N/m, k2 = 5 N/m, F = 3 N, and omega_f = pi rad/sec with all initial conditions equal to zero. Plot the two responses (x1, x2) from t = 0 to 50.

FYI, this system models two masses of unit dimensions (1kg) connected by a spring of stiffness k2, with one mass connected to ground through a spring of stiffness k1 being driven sinusoidally by a force of 3N.

Pendulum.m

h = 0.0001;t = [0:h:30];th = zeros(length(t),1);om = th;
th(1) = 0.3;
g = 9.80665;l = 1;
for i = 2: 1: length(t)    th(i) = th(i-1) + h*om(i-1);    om(i) = om(i-1) - h*g/l*sin(th(i-1));end
plot(t,th);hold on;

pend_driver.m

[T,OM] = ode23(@pend, [0 30], [0.3 0]);
plot(T, OM(:,1));hold on;

pend.m

function [ dydt ] = pend( t, y )

g = 9.80665;l = 1;
% dydt(1) = y(2);

% dydt(2) = -g/l*sin(y(1));

% dydt = dydt';
dydt = [y(2); -g/l*sin(y(1))];

end

Explanation / Answer

VBA to use Euler method

function uler(func,y0,delt,tf)

t=0:delt:tf;

y(1)=y0;

for i=1:length(t)-1

y(i+1)=y(i)+delt*(feval(func,t(i)));

end

plot(t,y);

xlabel('time')

ylabel(y')

Disp(y(end))

Considering the example

dy/dt = 14t-3

y(0)=5

delt=0.1

Now the equation function is as below in Macro editor

function z=eqn(t)

z=14*t-3;

end

Related Questions

Navigate

• Browse (All)
• Browse D
• Subjects

Integrity-first tutoring: explanations and feedback only — we do not complete graded work. Learn more.