Question

Using Euler method, ODE45 and Simulink solve the following ODE's.

a) y+c(dy/dt)=10sinwt

c=2, w=10, y(0)=0

Assume time span 0: 0.01: 10

b) dy/dt=t^2+yz

dz/dt=t+y^2z^2

y(0)=0.2

z(0)=-0.1

Assume time span 0: 0.1 :1

c)(d^2)y/dt^2+0.5(dy/dt)+siny=0

y(0)=1

dy/dt =0 when t=0

assume time span 0: 0.01 :10

Matlab code please

Answer 1

% PART (a):
% y + c dy/dt = 10 sin wt, y(0) = 0
% dy/dt = (10 sin wt - y) / c

clc,clear,close all

c = 2; w = 10; y0 = 0;
h = 0.01; % dt
t = 0 : h : 10; % time vector (span)
n = length(t) - 1; % number of sub-intervals

f = @(t,y) (10 * sin(w*t) - y) / c; % dy/dt as a function

% Euler's method
y(1) = y0; % set initial condition
for i = 1 : n
y(i+1) = y(i) + h * f(t(i),y(i));
end

% ode45 method
[~,y2] = ode45(f,t,y0);

% plot results
plot(t,y,'- r',t,y2,'- b')
xlabel('t'),ylabel('y(t)'),legend('euler','ode45')

% PART (b): dy/dt = t^2 + yz, y(0)=0.2, dz/dt = t + y^2z^2, z(0)=-0.1
clc,clear,close all

h = 0.1; t0 = 0; tn = 1;
t = t0 : h : tn;
y0 = 0.2; z0 = -0.1; % y(0) and z(0)
% write the odes as functions f(t,p) where p = [y,z] ie. y = p(1), z = p(2)
f1 = @(t,p) t^2 + p(1)*p(2); % dy/dt = t^2 + yz
f2 = @(t,p) t + p(1)^2 * p(2)^2; % dz/dt = t + y^2z^2

% Euler's method
y(1) = y0; z(1) = z0; % set the initial conditions
n = (tn - t0)/h;
for i = 1 : n
y(i+1) = y(i) + h * f1(t(i),[y(i) z(i)]);
z(i+1) = z(i) + h * f2(t(i),[y(i) z(i)]);
end

% ode45 method
f = @(t,p) [f1(t,p) ; f2(t,p)];
[~,P] = ode45(f,t,[y0 z0]);
y2 = P(:,1); % y solution
z2 = P(:,2); % z solution

% plot results
figure,plot(t,y,'-r',t,y2,'-b'),xlabel('t'),ylabel('y'),legend('euler','ode45')
figure,plot(t,z,'-r',t,z2,'-b'),xlabel('t'),ylabel('z'),legend('euler','ode45')

% PART (c):
% d^y/dt^2 + 0.5 dy/dt + sin y = 0, y(0) = 1, dy/dt(0) = 0
% y'' + 0.5y' + sin y = 0, y(0) = 1, y'(0) = 0
% let z = y'
% z' + 0.5z + sin y = 0, y(0) = 1, z(0) = y'(0) = 0
% z' = -0.5z - sin y
% solve the 2 odes simulteneously:
% y' = z, y(0) = 1
% z' = -0.5z - sin y, z(0) = 0

clc,clear,close all

t0 = 0; tn = 10;
h = 0.01;
t = t0 : h : tn; % time vector

% initial conditions y(0) and z(0)
y0 = 1; z0 = 0;   

% define the odes as functions f1(t,p) and f2(t,p) where p = [y,z]
f1 = @(t,p) p(2); % y' = z
f2 = @(t,p) -0.5*p(2) - sin(p(1)); % z' = -0.5z - sin y

% Euler's method
y(1) = y0; z(1) = z0; % set the initial conditions
n = length(t) - 1; % number of sub-intervals

for i = 1 : n
y(i+1) = y(i) + h * f1(t(i),[y(i) z(i)]);
z(i+1) = z(i) + h * f2(t(i),[y(i) z(i)]);
end

% ode45 method   
f = @(t,p) [ f1(t,p) ; f2(t,p)];   
[~,P] = ode45(f,t,[y0 z0]);
% extract solutions
y2 = P(:,1); % y solution

% plot results
figure,plot(t,y,'-r',t,y2,'-b')
xlabel('t'),ylabel('y(t)'),legend('euler','ode45')

------------------------------------------------------------------------------
COMMENT DOWN FOR ANY QUERY RELATED TO THIS ANSWER,

IF YOU'RE SATISFIED, GIVE A THUMBS UP
~yc~