% This program uses Euler's method to solve the initial value problem % y' = y(1-y) for the initial value (0,0.1). In this example the % differential equation is solved for t between 0 and 5 % function eulerexample2 tstart = 0; % Start time tend = 5; % End time dt = 0.1; % Delta t t = [0:dt:tend]; N = length(t); y0 = 0.1; % Initial condition yeu(1) = y0; % Apply the initial condition yan(1) = y0; yrk(1) = y0; % The general ODE: dy/dt = f(t,y) % Our example autonomous ODE: dy/dt = f(0,y) = y (1-y) % the small step size h = dt %Eulers method: y_n+1 = y_n + dt f(t_n, y_n) for i = 2:N yeu(i) = yeu(i-1) + fct(yeu(i-1))*dt; end; %2nd order RK method: y_n+1 = y_n + dt f(t_n + dt/2, y_n + dt/2 f(t_n, y_n)) %Our example: y_n+1 = y_n + dt f(y_n + dt/2 f(y_n) ) since we have no explicity t dependence for i = 2:N K1 = dt * fct(yrk(i-1)); K2 = dt* fct (yrk(i-1) + 1/2 * K1); yrk(i) = yrk(i-1) + K2; end; % Solved numerically using matlab options = odeset('RelTol',1e-6,'Abstol',1e-6); sol = ode45(@rhs,[tstart tend],y0,options); % function name, time interval, initial conditions, options, e.g. parameters or other arguments yod = deval(sol,t); % interpolates the solution at t = 0 0.1 0.2 0.3 .... % Analytical solution: yan = 1./(1 + 9*exp(-t)); % Error Eeu = abs(yeu-yan); Erk = abs(yrk-yan); Eod = abs(yod-yan); figure(3); h=plot(t,yan,'k',t,yeu,'b',t,yrk,'r',t,yod,'c'); set(h,'linewidth',2); set(gca,'fontsize',20); xlabel('time'); ylabel('y'); legend('ana','eul','rk2','ode45'); xlim([0 5]); grid; figure(4); h=plot(t,Eeu,t,Erk,'r',t,Eod); set(h,'linewidth',2); set(gca,'fontsize',20); xlabel('time'); ylabel('error'); legend('euler','rk2','ode45'); xlim([0 5]); grid; figure(5); h=plot(t,Erk,t,Eod,'r'); set(h,'linewidth',2); set(gca,'fontsize',20); xlabel('time'); ylabel('error'); xlim([0 5]); grid; figure(5); h=plot(t,Eod,'r'); set(h,'linewidth',2); set(gca,'fontsize',20); xlabel('time'); ylabel('error'); xlim([0 5]); grid; function dy = rhs(t,y) dy = y*(1-y); function f = fct(y) f = y*(1-y);