Question

Assignment 9.3: ODE Ill or pendulum Ill Solve the equation for the small amplitude pendulum without driving forces but with damping: dt2 Use your RK4 routine. Plot the result.

using Matlab

Answer 1

Int we have to np+ itinto toso i^t or des-egm Let Initial comditons dtltoo dar So that diCo) a d 2- d toLemalab fn naw Here are will aing K4 The input vaxiables are- iuiti al t tin tend Final cin initial da din dt step size fmcti on 2 (o parametes

clear all close all %all parameter values al0 0; a20-5; tin-0; t end - 10 ; %initial and final time h 0.005; Binitial quess tep size all parameter values =1 ;q =.5; t rk,al_rk, a2_rk]= RK4_sol(t_in,t_end,a10,a20,h,g,1,q); Plotting data figure(1) plot (t_rk,al_rk) title'time vs. Displacement plot') xlabeltime') ylabeldisplacement') figure (2) plot (t rk, a2_rk) title'time vs. Velocity plot') xlabel'time' ylabelvelocity') %Function for RK4 solution function [t_rk,al_rk,a2_rk]= RK4_sol (t_in,t_end, al0, a20, h, g,1,q) function forRK4 equation solution fl-e (t,al, a2 ,g,l,q) a2; f2- (t,al, a2 ,g,1,q) (-g/1) *al-q*a2; n (t end-t_in)/h; al_rk1) a10; a2_rk(1)=a20; t_rk(1)=t_in; %number of steps %Runge Kutta 4 iterations for i-1:n kl-h*f1 (t rk (i),al_rk(i),a2_rk(i),g,1,q); 11-h*f2(t_rk (i), al_rk(i),a2_rk(i),g,1,q); k2-h*f1 (t_rk(i) +h/2,al_rk (i)+ (1/2) *kl,a2_rk(i)+1/2) l1,g,1,q); 12 h*f2 (t_rk (i)+h/2,al_rk(i)+(1/2)kl,a2_rk(i)+(1/2) 11,g,l,q); k3-h*f1 (t_rk(i) +h/2,a1_rk(i)+(1/2)*k2,a2_rk(i)+ (1/2 ) 12,g,1,q); 13-h*f2 (t_rk (i) +h/2,al_rk(i)+(1/2 ) * k2,a2_rk(i)+ (1/2 ) 12,g,1,q); k4 h*f1(t_rk(i) +h, al_rk (i) +k3, a2_rk (i) +13,g,1,q); 14-h*f2(t_rk (i) +h, al_rk(i) +k3, a2_rk (i) +13,g,1,q); t_rk(i+1)t_in+i *h; al_rki+1=double(al_rk(i)+1/6) * (kl+2 *k2+2*k3+k4)); a2_rk(i+1 double(a2_rk(i)+( 1/6) * ( 11+2 *12+2*13+14) ); end end 웅웅응용웅용용음용용음용융용용용용용용 End 융응음용용용용용용용용용8용용용8용용 of Code 1

time vs. Displacement plot 15 0.5 0.5 10 1.5 6 5 time time vs. Velocity plot 3. 10 6 5 time 2

clear all
close all
%all parameter values
a10=0; a20=5;     %initial guess
t_in=0; t_end=10; %initial and final time
h=0.005;          %step size
g=9.8; l=1;q=.5; %all parameter values
%all solution plot
[t_rk,a1_rk,a2_rk]= RK4_sol(t_in,t_end,a10,a20,h,g,l,q);
%Plotting data
figure(1)
plot(t_rk,a1_rk)
title('time vs. Displacement plot')
xlabel('time')
ylabel('displacement')

figure(2)
plot(t_rk,a2_rk)
title('time vs. Velocity plot')
xlabel('time')
ylabel('velocity')
%Function for RK4 solution
function [t_rk,a1_rk,a2_rk]= RK4_sol(t_in,t_end,a10,a20,h,g,l,q)
    %function forRK4 equation solution
    f1=@(t,a1,a2,g,l,q) a2;
    f2=@(t,a1,a2,g,l,q) (-g/l)*a1-q*a2;
  
        n=(t_end-t_in)/h; %number of steps
        a1_rk(1)=a10; a2_rk(1)=a20; t_rk(1)=t_in;
    %Runge Kutta 4 iterations
        for i=1:n
            k1=h*f1(t_rk(i),a1_rk(i),a2_rk(i),g,l,q);
            l1=h*f2(t_rk(i),a1_rk(i),a2_rk(i),g,l,q);
            k2=h*f1(t_rk(i)+h/2,a1_rk(i)+(1/2)*k1,a2_rk(i)+(1/2)*l1,g,l,q);
            l2=h*f2(t_rk(i)+h/2,a1_rk(i)+(1/2)*k1,a2_rk(i)+(1/2)*l1,g,l,q);
            k3=h*f1(t_rk(i)+h/2,a1_rk(i)+(1/2)*k2,a2_rk(i)+(1/2)*l2,g,l,q);
            l3=h*f2(t_rk(i)+h/2,a1_rk(i)+(1/2)*k2,a2_rk(i)+(1/2)*l2,g,l,q);
            k4=h*f1(t_rk(i)+h,a1_rk(i)+k3,a2_rk(i)+l3,g,l,q);
            l4=h*f2(t_rk(i)+h,a1_rk(i)+k3,a2_rk(i)+l3,g,l,q);
            t_rk(i+1)=t_in+i*h;
            a1_rk(i+1)=double(a1_rk(i)+(1/6)*(k1+2*k2+2*k3+k4));
            a2_rk(i+1)=double(a2_rk(i)+(1/6)*(l1+2*l2+2*l3+l4));
        end
end

%%%%%%%%%%%%%%%%%%%% End of Code %%%%%%%%%%%%%%%%%%%