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

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

Answer 1

[Roger Stafford]

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.

[Laura]

1) You meant put x=xin under for dt loop?

4) I want to use RK4 and Euler to see what the errors are.

[Roger Stafford]

1) Yes, the outer loop with changing dt. The trouble with that is that you will lose all the information about your data for dt - .0001 . You need to store it somewhere if you want to do calculation with it later.