clear; format long; format compact % Recall from Maple % numbers:={gamma=1, omega=sqrt(5/4), Theta=Pi/4, Omega=2}; gamma=1; omega=sqrt(5/4); Theta=pi/4; Omega=2; T=15; % Time to stop integrating % Recall ODE from Maple script: % diff(theta(t),t,t) + gamma*diff(theta(t),t) + omega^2*sin(theta(t))=0 % Create a function (handle) giving the ODE to be solved: dy_dt = @(t,y) [y(2); -gamma*y(2) - omega^2 * sin(y(1)) ]; % Now solve the ODE using Matlab's built in function (same as Maple): [t,Y] = ode45(dy_dt, [0 T], [Theta Omega]); n_steps = size(t,1) figure(1); clf plot(t,Y(:,1),'o--r', t,Y(:,2),'o--b'); title(['Numerical solution']); figure(2); clf plot(Y(:,1), Y(:,2), 'o--r'); title(['Phase trajectory']); %return % Now let us try to do this ourselves using 100 steps of the Euler method: N_steps=100; % Number of steps (play around with this) dt = T/N_steps; % Time step size t_E = zeros(N_steps+1,1); % Initialize memory for t's Y_E = zeros(N_steps+1,2); % Initialize memory for Y's t_E(1)=0; Y_E(1,1)=Theta; Y_E(1,2)=Omega; % Initial conditions for step=2:N_steps+1 t_E(step) = t_E(step) + dt; % Time Y_E(step,:) = Y_E(step-1,:) + dt * dy_dt(t_E(step-1), Y_E(step-1,:))'; end figure(1); hold on; plot(t_E,Y_E(:,1),'sg', t_E,Y_E(:,2),'sk'); legend(['rk45: theta', 'rk45: omega', 'Euler: theta', 'Euler: omega']); figure(2); hold on; plot(Y_E(:,1), Y_E(:,2), 's--b'); legend(['rk45', 'Euler']);