%==================================================================== % Initial value problem for a single 1st-order ODE % Call three methods: Euler, Predictor-Corrector, Runge-Kutta 4th order %==================================================================== clear; clc; %open a file if needed % fileID = fopen('results.txt','w'); % initial conditions xi = 1.00; ti = 0.00; % "time" step and max time dt = 0.09; tmax = 10.0; %select function dx = @f; %preparing to plot i = 1; t (i)= ti; x1(i) = xi; x2(i) = xi; x3(i) = xi; % print the header and initial conditions % fprintf(' Solving 1st order single ODE \n') % % fprintf(' t x(t) \n') % formatSpec1 = '%10.5f %10.5f \n'; % fprintf(formatSpec1,ti,xi) %fprintf(fileID,formatSpec1,ti,xi) xi1 = xi; xi3 = xi; while(ti < tmax) tf = ti + dt; i = i+1; [xf1] = euler1(dx,ti,xi1,tf); x1(i) = xf1; xi1 = xf1; [xf3] = RK4D11(dx,ti,xi3,tf); x3(i) = xf3; xi3 = xf3; t(i) = tf; ti = tf; % if(key==2) % [xf] = euler2(dx,ti,xi,tf); % end % if(key==3) % [xf] = RK4D11(dx,ti,xi,tf); % end % xtrue = exp(-1.0*tf); % fprintf(formatSpec1,tf,xf); % fprintf(fileID,'%10.5f %10.5f %10.5f\n',tf,xf,xtrue); % to plot end %while plot(t,x1,'b',t,x3,'r') function f=f(t,x) %======== % function for dx/dt = f(t,x) %======== % f = (-1.0)*x+t*0.0; f = (-1.0)*x*cos(t)+0.5*sin(2*x); % f = t*x*x; end function [xf] = euler1(dx,ti,xi,tf) %================================================== % Solution of a single 1st order ODE dx/dt=f(x,t) % Method: Simple Euler % Alex G. October 2020 %-------------------------------------------------- % input ... % dx(t,x)- function dx/dt (supplied by a user) % ti - initial time % xi - initial position % tf - time for a solution % output ... % xf - solution at point tf %================================================== xf = xi + dx(ti,xi)*(tf-ti); end function [xf] = euler2(dx,ti,xi,tf) %================================================== % Solution of a single 1st order ODE dx/dt=f(x,t) % Method: Euler: predictor-corrector % Alex G. October 2020 %-------------------------------------------------- % input ... % dx(t,x)- function dx/dt (supplied by a user) % ti - initial time % xi - initial position % tf - time for a solution % output ... % xf - solution at point tf %================================================== xf = xi + dx(ti,xi)*(tf-ti); xf = xi + (dx(ti,xi)+dx(tf,xf))*0.5*(tf-ti); end function [xf] = RK4D11(dx,ti,xi,tf) %================================================== % Solution of a single 1st order ODE dx/dt=f(x,t) % Method: Runge-Kutta 4th order % Alex G. October 2020 %-------------------------------------------------- % input ... % dx(t,x)- function dx/dt (supplied by a user) % ti - initial time % xi - initial position % tf - time for a solution % output ... % xf - solution at point tf %================================================== h = tf-ti; k1 = h*dx(ti,xi); k2 = h*dx(ti+h/2.0,xi+k1/2.0); k3 = h*dx(ti+h/2.0,xi+k2/2.0); k4 = h*dx(ti+h,xi+k3); xf = xi + (k1 + 2.0*(k2+k3) + k4)/6.0; end