Euler's Method for stiff ODE

Having problem with my code, I'm trying to make my graph (the middle) look like the bottom graph. I feel like it is something with the initial inputs, I pretty sure my for loop is correct. Thanks for the help I appreciate it!

if true

    % %%Explicit Euler
      h = 5;
      x = [0:h:20];
      yb = zeros(size(x));
      yb(1) = 1;
      yb2(1) = 1;
      mu = 1;

       for n = 1:length(x)-1
         yb(n+1) = yb(n) + h*yb2(n);
         yb2(n+1) = yb2(n) + h*(mu*(1-yb(n)^2)*yb2(n) - yb(n));
      end

      %%Plot
      figure(1); clf(1);
      plot(x,yb,'k');
      hold on
      plot(x,yb2,'or');
      title(['Van der Pol Equation for \mu = ',num2str(mu), 'with Euler']);
      xlabel('Time'); 
      ylabel('y(t)');
      legend('y_{1}', 'y_{2}','location','northeast')

end

[Torsten]

Honestly: How do you want to reproduce the curvature of the bottom graph with a step size of h=5 ???

[jake thompson]

That's what I figured also but everytime I change my step size every thing gets a bit strange.

[Torsten]

Could you show the plot with stepsize h=0.01, say ?

[Rena Berman]

(Answers Dev) Restored edit

Answer 1

[Jim Riggs]

Torsten is right. There is no way you will get a solution with a large step size (h=5). I checked out your code and it works fine for mu=1 with a step size of h=0.1

I would recommend the following two changes;

1) add the command

 clear x yb yb2

at the start of the script.

2) I would also recommending pre-allocating the yb2 vector. you can do this by changing

 yb2(1) = 1

to

 yb2 = yb.

this will make yb2 and yb the same size, and they will both have a 1 in the first position.

Now, as you increase the value of mu, the system becomes more stiff. For example, using a step size of h=0.1 I can get a solution for mu =2 and mu=3, but at mu=4 the function blows up after about 14 seconds. (it is numerically unstable). Now you have to reduce the step size to improve the calculation stability. In order to get a solution with mu=1000, and to execute for 3000 seconds without blowing up, the step size must be extremely small (10^-6 or even smaller). My problem is that this is taxing the capabilities of my computer and I cannot find the solution.

A major part of the problem is the use of the Euler integration method. This method is not well suited for stiff systems (perhaps that is what this lesson is about). For this problem, I would switch to a Runge-Kutta integrator and try my favorite [5,6] Runge Kutta method. This would allow the use of a much larger step size (probably 3 orders of magnitude) compared to what is required for Euler integration.

[jake thompson]

I edited the first few lines of my code and solved my problem!

if true
   m = 500; %  number of steps
   t_span = [0 20]; %time span
   h = t_span(end)/(m-1);
   x = [0:h:t_span(end)];
end

[Jim Riggs]

This will work for mu=1, since it results in a step size of 0.04. Bit it will not work for significantly larger values of mu. As mu increases, you still need to decrease the step size. You will have to experiment to find acceptable values. And when you increase the timespan, the step size gets bigger. This is also not desirable.