Plot Euler's method on Matlab

Viewed 597

I made the code for euler's method in matlab and now I have to plot the approximation and the exact result. The problem is that I don't know how to introduce the analytical solution and plot it. I made this but it doesn't work.

function [t,w] = euler (f,y0,a,b,h) %func=f(x,y)=dy/dx               
                                            %a and b the interval ends 
                                            %h=distance between partitions
                                             %y0= initial value
                s=input('Give an equation in x: ','s'); %the solution of the ODE
                n=(b-a)/h;
                t=zeros(1,n+1);
                w=zeros(1,n+1);
                t(1)=a;
                w(1)=y0;

              for i=1:n
                    w(i+1)=w(i)+h*f(t(i),w(i));
                    t(i+1)=a+h*i;

              end
          w
          x=a:h:b;  
          y=s;              
          plot(x,y,'k')
          hold on             
          plot(t,w,'r')      
           hold off

endfunction

The analytical solution is y=-x-1+2*exp(x) and the ODE is y'=x+y.

1 Answers

the explicit euler method is not unconditionally stable, so make sure, that you pick a sufficiently small timesteps for your tests. I barely changed anything in your code, but I needed to make it slightly readable. If the problem persists, pls upload a screenshot of the plot.

function [x,y] = euler(y0,a,b,h)        

      % Left out f for algorithm test!            
      % a and b the interval ends
      % h=distance between partitions ~ 
      % y0= initial value

      s=@(x) 2*exp(x) - (x+1);
      f=@(x,y) x + y;
      %s=input('Give an equation in x: ','s'); %the solution of the ODE
      x = a:h:b;
      n=length(x);
      y=zeros(n,1);
      y(1)=y0;

      for i=1:n-1
          y(i+1)=y(i)+h*f(x(i),y(i));
      end

      plot(x,s(x),'k')
      hold on             
      plot(x,y,'r')      
      hold off
endfunction

I'd recommend you sticking to the least amount of variables as possible. Try not mixing the time and space domain, often functions depend on both and it is then confusing to have a y(x) with timesteps.

Img_Plot

It should work out now. Notice that y0 has to match a point of your solution function on a given value for a.

s(a) = y0

Cheers,

Pablo

Related