Question

Runge-Kutta method R-K method is given by the following algorithm. Yo = y(xo) = given. k1-f(xy) k4-f(xi +h,yi + k3) 6 For i = 0, 1, 2, , n, where h = (b-a)/n. Consider the same IVP given in problem 2 and answer the following a) Write a MATLAB script file to find y(2) using h = 0.1 and call the file odeRK 19.m b) Generate the following table now using both ode Euler and odeRK19 only for h -0.01. yi (EM) yi (RK) Exact yi Error-EM (%) Error-RK(%) Xi 0.0 0.5 3.000000 4.080575 0.00 0.02 3.000000 3.000000 4.072295 0.00 2.0 c) Use the subplot command to produce plots for exact solution, EM and RK method for h 0.1 and h-0.01 in two different figures. (That means in one graph they have plots for the exact solution, EM and RK for h -0.1 and in another graph for h- 0.01) d) What can you say about the convergence rate of EK and RK method according to part (c)? dy dr =-1.2y + 7e-0.3x from x = 0 to x = 2.0 with the initial condition y-3 atx-0

Answer 1

SMatlab code for Euler's forward and Rk4 method clear all close all %A11 h values hh 0.1,0.01 for it-1:2 Program for Euler %function for Euler equation solution f-L(x,y) -1.2*y+7.*exp(-0.3*x) %Euler steps for h and x-end Initial values y0-3 ; %step size h-hh (it) %x end value xend-2; % Euler steps y_resultl (1)-y0; x_resultl (1)-x0; %Loop for Euler Steps for i-1: length (xn)-1 x_resultl (i+1)x_resultl (i) +h; y resultl(i +1)-y_resultl(i) +h*double(f (x resultl(i),y resultl (i))); end printing the result for Euler fprintf( "\tThe approximate solution using Euler method at x=%d for h-% 0 . 2f is %f\n" , xe nd , h , y-result1(end )) Program for RK4 %function for RK 4 eguation solution f-ê (x,y) -1.2*y+7.*exp (-0.3*x); RK4 steps for h and x end %Initial values x0 0 y0-3 ; end value xend-2 xn x0:h: xend; % RK4 steps lt2 ( 1 ) =yo; y resu xresult2(1)#x0; - %Runge Kutta 4 iterations n-(xend-x0)/h; for i-l:n

k0sh* f (x-result2( ǐ ) , y-result2( i)) ; k1-h*f(x-result2(1) + (1/2) *h, y-re sult 2 (i) +(1/2) *ko) ; k2-h * f(x-result2(ǐjt (1/2) *h,y-result2(i) + (1/2) *k1); k3-h * f (x-result2 ( İ )th, y-result 2 ( i)-k2 ) ; x result2 (i+1)x result2 (i) +h; y-result 2 ( İ +1)-double (y result2(i)+ (1/6)* (k0+2*kl+2 k2+k3)); end tprinting the result for RK4 fprintf("1tThe approximate solution using RK4 method at x-2d for h=% 0 . 2f is %f\n.xe nd , h, y-result2 ( end )) %Finding exact solution using dsolve syms y(x) eqn-diff (y,x)1.2*y+7.*exp(-0.3*x); cond y (0)3; y_result3(x)-dsolve eqn, cond); printing exact solution fprintf( "\tThe exact solution using dsolve at x-%d is %f n',xend,y result3 (xend) ) fprintf ( The expression for Exact solution\n') disp(y_result3 (x)) figure (it) subplot(2,1,1) hold on plot (x_resultl,double(y result3(x_result1))) plot (x_result1,y_resultl) legend ( Exact solution', 'Euler Solution') title( sprint f (Euler and Exact solution for h xlabel('x) ylabel'y') 2.3f',h)) subplot (2,1,2) hold on plot (x_resultl,double(y result3(x resultl))) plot (x_result2,y result2) legend( Exact solution, 'RK4 Solution') title(sprintf("RK4 and Exact solution for h-%2.3f',h)) xlabel(x) ylabel'y) if it1 Printing the result for table fprint f ( 'Table for h-o. 1 \n") fprintf n\tx,\ty (EM),Nty (RK),\tExact (Y),\tError-EM,tError- for ii-l:length(y resultl) ext_sl-y_result3 (x_resultl(ii)); err1-double( (abs (y resultl(ii)-ext sl)/ext sl) *100):

err2-double( (abs (y result2 (ii)-ext sl)/ext_s1) 100); x result1(ii),y_resultl(ii),y result2 (ii),ext_sl,errl,err2) end end 용용号号%용号응용용号응용용융융용 End of Code 융융욤욤융융욤욤융융욤욤응응욤욤 The approximate solution using Euler method at x=2 for h=0.10 is 3.879902 The approximate solution using RK4 method at x=2 for h=0.10 is 3.835103 The exact solution using dsolve at x=2 is 3.835105 The expression for Exact solutiorn (70*exp(-(6*x)/5) *exp( (9*x)/10))/9 - (43 texp (- (6*x)/5))/9 Table for h=0.1 x, y(EM), y (RK),Exact(Y), Error-EM, Error-RK, 0.00, 3.00000, 3.00000, 3.00000, 0.00, 0.00000, 0.10, 3.34000, 3.31040, 3.31040, 0.89, 0.00003, 0.20, 3.61851, 3.56650, 3.56650, 1.46, 0.00004, 0.30, 3.84353, 3.77501, 3.77501, 1.81, 0.00005, 0.40, 4.02205, 3.94186, 3.94186, 2.03, 0.00006, 0.50, 4.16025, 4.07229, 4.07230, 2.16, 0.00006, 0.60, 4.26352, 4.17095, 4.17095, 2.22, 0.00007, 0.70, 4.33658, 4.24192, 4.24193, 2.23, 0.00007, 0.80, 4. 38360, 4.28884, 4.28884, 2.21, 0.00007, 0.90, 4.40821, 4.31488, 4.31488, 2.16, 0.00006, I.00, 4.41359, 4. 32288, 4.32288, 2.10, 0.00006, 1.10, 4.40253, 4.31531, 4.31532, 2.02, 0.00006, 1.20, 4.37748, 4.29438, 4.29438, 1.93, 0.00006, 1.30, 4.34055, 4.26201, 4.26201, 1.84, 0.00005, 1.40, 4.29363, 4. 21991, 4.21991, 1. 75, 0.00005, 1.50, 4.23832, 4.16957, 4.16957, 1.65, 0.00005, I.60, 4.17606, 4.11230, 4.11230, 1.55, 0.00005, 1.70, 4.10808, 4.04927, 4.04927, 1.45, 0.00004, I.80, 4.03546, 3.98149, 3.98149, 1. 36, 0.00004, 1.90, 3.95913, 3.90984, 3.90984, 1.26, 0.00004, 2.00, 3.87990, 3.83510, 3.83510, 1. 17, 0.00004, The approximate solution using Euler method at x=2 for h=0.01 is 3.839567 The approximate solution using RK4 method at x-2 for h-0.01 is 3.835105 The exact solution using dsolve at x-2 is 3.835105 The expression for Exact solution (70*exp(-(6*x)/5) *exp( (9*x)/10)/9 - (43 *exp(- (6*x)/5))/9

Euler and Exact solution for h-0.100 4.5 -Exact solution Euler Solution 3.5 0.5 1.5 25 RK4 and Exact solution for h-0.100 Exact solution RK4 Solution 3.5 0.5 1.5 25 Euler and Exact solution for h:0.010 4.5 Exact solution -Euler Solution 3.5 0.5 1.5 25 RK4 and Exact solution for h-0.010 4.5 ㅡ Exact solution RK4 Solution 3.5 0.5 1.5 25

%%Matlab code for Euler's forward and Rk4 method
clear all
close all

%---------------------%
%All h values
hh=[0.1,0.01];
for it=1:2
    %Program for Euler
    %function for Euler equation solution
    f=@(x,y) -1.2*y+7.*exp(-0.3*x);
    %Euler steps for h and x_end
    %Initial values
    x0=0;
    y0=3;
    %step size
    h=hh(it);
    %x end value
    xend=2;
    xn=x0:h:xend;
    % Euler steps
    y_result1(1)=y0;
    x_result1(1)=x0;
    %Loop for Euler Steps
    for i=1:length(xn)-1
        x_result1(i+1)= x_result1(i)+h;
        y_result1(i+1)=y_result1(i)+h*double(f(x_result1(i),y_result1(i)));
    end
    %printing the result for Euler
    fprintf('\tThe approximate solution using Euler method at x=%d for h=%0.2f is %f\n',xend,h,y_result1(end))
  
    %---------------------%
    %Program for RK4
    %function for RK4 equation solution
    f=@(x,y) -1.2*y+7.*exp(-0.3*x);
    %RK4 steps for h and x_end
    %Initial values
    x0=0;
    y0=3;
    %x end value
    xend=2;
    xn=x0:h:xend;
    % RK4 steps
    y_result2(1)=y0;
    x_result2(1)=x0;

        %Runge Kutta 4 iterations
        n=(xend-x0)/h;
        for i=1:n

            k0=h*f(x_result2(i),y_result2(i));
            k1=h*f(x_result2(i)+(1/2)*h,y_result2(i)+(1/2)*k0);
            k2=h*f(x_result2(i)+(1/2)*h,y_result2(i)+(1/2)*k1);
            k3=h*f(x_result2(i)+h,y_result2(i)+k2);
            x_result2(i+1)=x_result2(i)+h;
            y_result2(i+1)=double(y_result2(i)+(1/6)*(k0+2*k1+2*k2+k3));

        end
    %printing the result for RK4
    fprintf('\tThe approximate solution using RK4 method at x=%d for h=%0.2f is %f\n',xend,h,y_result2(end))

    %---------------------%
    %Finding exact solution using dsolve
    syms y(x)
    eqn=diff(y,x)==-1.2*y+7.*exp(-0.3*x);
    cond=y(0)==3;
    y_result3(x)=dsolve(eqn,cond);
    %printing exact solution
    fprintf('\tThe exact solution using dsolve at x=%d is %f\n',xend,y_result3(xend))
    fprintf('The expression for Exact solution\n')
    disp(y_result3(x))

    figure(it)
    subplot(2,1,1)
    hold on
    plot(x_result1,double(y_result3(x_result1)))
    plot(x_result1,y_result1)
    legend('Exact solution','Euler Solution')
  
    title(sprintf('Euler and Exact solution for h=%2.3f',h))
    xlabel('x')
    ylabel('y')
  
    subplot(2,1,2)
    hold on
    plot(x_result1,double(y_result3(x_result1)))
    plot(x_result2,y_result2)
    legend('Exact solution','RK4 Solution')

    title(sprintf('RK4 and Exact solution for h=%2.3f',h))
    xlabel('x')
    ylabel('y')
  
    if it==1
        %Printing the result for table
        fprintf('Table for h=0.1\n')
        fprintf('\n\tx,\ty(EM),\ty(RK),\tExact(Y),\tError-EM,\tError-RK,\n')
        for ii=1:length(y_result1)
            ext_sl=y_result3(x_result1(ii));
            err1=double((abs(y_result1(ii)-ext_sl)/ext_sl)*100);
            err2=double((abs(y_result2(ii)-ext_sl)/ext_sl)*100);
            fprintf('\t%2.2f,\t%2.5f,\t%2.5f,\t%2.5f,\t%2.2f,\t%2.5f,\t\n',...
            x_result1(ii),y_result1(ii),y_result2(ii),ext_sl,err1,err2)
        end
    end

end

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