MATLAB: Calculate the error for many time steps in the same file ( Only work for one time step )

errorstore data

I want to calculate the error for my data at different time step (dt) in one file. If I tried one time step, dt, then it works. However, it does not work for many time steps at once. For example, if I have dt=0.0001:0.0001:0.0001, it will give me the right answer. However, if I use dt=0.0001:0.0001:0.0002, then the first data for dt=0.0001 is correct, but the data for dt=0.0002 is incorrect.
 % Solve dx/dt=H(x)+1
 % H(x)=0 when 0<=x<=0.1
 % H(x) = 1, when x>0.1
 % Time step
 First file 
 function x=heaviside(x,xo)
 if x<xo
    x=1;
 else x>xo
    x=2;
 end
 end
Second file
 j=0;
xin=0; % Initial position at t=0
xo=0.1;
x=xin;
for dt=0.0001:0.0001:0.0002
j=j+1;
tf=0.2; % final time
t=[0:dt:tf];
n=length(t);
 for i=1:n-1
    % RK4 method
   k1=heaviside(x,xo);
   k2=heaviside(x+(0.5*dt*k1)',xo);
   k3=heaviside(x+(0.5*dt*k2)',xo);
   k4=heaviside(x+(dt*k3)',xo);
   x(i+1,:) = x(i,:) + ((k1+2*k2+2*k3+k4)/6)'*dt;
 end
 end
Error(j)=abs(((x(n,:)-0)-(0.3))/0.3); % Calculate the error 
clear n j dt t x xo
where, x(n,:) is the final point of x vector, 0 is initial, 0.3 is exact solution. Please helps. Thanks

Best Answer

Questionable aspects of this code:
1) In your outer for-loop you don't reset x to an initial value x = xin again for a rerun with different dt. This would show up with dt = 0.0002 .
2) You appear to be calling on the 'heavyside' function with x a multiple-element vector, and yet it is not designed to handle such a vector since it applies the if-else-end coding to it. 'if' applied to a logical vector will be true only if all elements of the vector are true. I don't think this is what you want.
3) On the other hand, in the Runge Kutta method I don't think 'heavyside' should be called upon with a vector x. It should be a scalar quantity, namely x(i).
4) The RK4 Runga Kutta method is designed to do fourth order approximation to the integrand function involved, but in this case your function has all its derivatives discontinuous at the point xo, so this might result in greater errors than you expect from the method. It is not designed to handle such functions.
Related Question