Question

Numerical Methods for Differential Equations - Please post full correct solution!!! - need to use MATLAB

3. (a) Write Matlab functions to integrate the initial value problem y = f(x,y), y(a) = yo, on an interval [a, b] using: • Euler's method • Modified Euler • Improved Euler • Runge Kutta 4 It is suggested that you implement, for example, Improved Euler as [x, y) = eulerimp('f', a, yo, b, stepsize), where (2,y) = (In, Yn) is the computed solution. (b) Use the four Matlab function you wrote in (a) with stepsize h = 1/40 and 1/80 to solve the following ODE: 1 V = -v*/2, y(0)=1 with Yexact (2) = (2) Calculate the absolute errors (i.e. Yezact (2) - Ycomputed(x)) at 2 = 1 to see if they are roughly reduced by a half for Euler's method, by a quarter for Modified Euler and Improved Euler, and by 1/16 for RK4, as the stepsize h is halved from 1/40 to 1/80. Discuss the performance of the methods by comparing the computation time used by each method to calculate y(1) with stepsize h = 1/40 and 1/80 (use the Matlab functions tic and toc to do the timing).

Answer 1

We have developed a MATLAB code for the given problem.

-------------------------MATLAB Code-------------------------

f=@(x,y) -0.5*y.^3; % f(x,y)
fe=@(x) 1./(sqrt(1+x)); % exact function
a=0;
y0=1;
b=1;
stepsize=[1/40 1/80];
yexact=fe(b);

for i=1:2
h=stepsize(i);
tic;[x1,y_eul]=euler(f,a,y0,b,h);T1=toc;
fprintf('\n\nEuler method: stepsize=%g\n',h);
fprintf('---------------------------------------------\n');
fprintf('y(%g)=%.8f\t error=%.10f\t Time=%g sec\n',b,y_eul,abs(y_eul-yexact),T1);
tic;[x2,y_mod]=eulermodified(f,a,y0,b,h);T2=toc;
fprintf('\n\nModified Euler method: stepsize=%g\n',h);
fprintf('---------------------------------------------\n');
fprintf('y(%g)=%.8f\t error=%.10f\t Time=%g sec\n',b,y_mod,abs(y_mod-yexact),T2);
tic;[x3,y_imp]=eulerimproved(f,a,y0,b,h);T3=toc;
fprintf('\n\nImproved Euler method: stepsize=%g\n',h);
fprintf('---------------------------------------------\n');
fprintf('y(%g)=%.8f\t error=%.10f\t Time=%g sec\n',b,y_imp,abs(y_imp-yexact),T3);
tic;[x4,y_rk4]=RK4(f,a,y0,b,h);T4=toc;
fprintf('\n\nRK4 method: stepsize=%g\n',h);
fprintf('---------------------------------------------\n');
fprintf('y(%g)=%.8f\t error=%.10f\t Time=%g sec\n',b,y_rk4,abs(y_rk4-yexact),T4);
end

function [x,y]=euler(f,a,y0,b,stepsize)
  
h=stepsize;
xn=a:h:b;
yn=0*xn; % initializing the array;
yn(1)=y0;
  
for i=1:numel(xn)-1
yn(i+1)=yn(i)+h*f(xn(i),yn(i));
end
x=xn(end);
y=yn(end);
end

function [x,y]=eulermodified(f,a,y0,b,stepsize)
  
h=stepsize;
xn=a:h:b;
yn=0*xn; % initializing the array;
yn(1)=y0;
  
for i=1:numel(xn)-1
yp=yn(i)+h*f(xn(i),yn(i));
yn(i+1)=yn(i)+h*f(xn(i)+0.5*h,yp);
end
x=xn(end);
y=yn(end);
end

function [x,y]=eulerimproved(f,a,y0,b,stepsize)
  
h=stepsize;
xn=a:h:b;
yn=0*xn; % initializing the array;
yn(1)=y0;
  
for i=1:numel(xn)-1
yp1=f( xn(i),yn(i) );
yp2=f( xn(i+1),yn(i)+h*yp1 );
yn(i+1)=yn(i)+0.5*h*( yp1+yp2 );
end
x=xn(end);
y=yn(end);
end

function [x,y]=RK4(f,a,y0,b,stepsize)
  
h=stepsize;
xn=a:h:b;
yn=0*xn; % initializing the array;
yn(1)=y0;
  
for i=1:numel(xn)-1
k1 = f(xn(i),yn(i));
k2 = f(xn(i)+0.5*h,yn(i)+0.5*h*k1);
k3 = f((xn(i)+0.5*h),(yn(i)+0.5*h*k2));
k4 = f((xn(i)+h),(yn(i)+k3*h));
  
yn(i+1) = yn(i) + (1/6)*(k1+2*k2+2*k3+k4)*h;
end
x=xn(end);
y=yn(end);
end

-------------------------MATLAB Output-------------------------

Euler method: stepsize=0.025 y (1)=0.70478126 error=0.0023255216 Time=0.0021672 sec Modified Euler method: stepsize=0.025 - - - - - - - - - - - - - y (1)=0.70950296 error=0.0023961751 Time=0.0014727 sec Improved Euler method: stepsize=0.025 y (1)=0.70712073 error=0.0000139537 Time=0.0020431 sec RK4 method: stepsize=0.025 y (1)=0.70710678 error=0.0000000001 Time=0.003726 sec

Euler method: stepsize=0.0125 -- - - y (1)=0.70595109 error=0.0011556953 Time=0.0006278 sec Modified Euler method: stepsize=0.0125 y (1)=0.70827994 error=0.0011731629 Time=0.0010728 sec Improved Euler method: stepsize=0.0125 - - - - - - - - --- y(1)=0.70711025 error=0.0000034707 Time=0.0009054 sec RK4 method: stepsize=0.0125 - - - - - - - - - y (1)=0.70710678 error=0.0000000000 Time=0.0019064 sec