How can I fix this?

   %Euler method
  %Forward, Backward, Modified, Predictor-Corrector, Analytical Solution, and Absolute Error in a table using %fprintf function.
  %dy/dx=2xy from x=1 to x=1.5. y(1)=1 stepsize is 0.05;

dx=0.05;
nfinal=0.5;
nsteps=(nfinal/dx);
xi=zeros(1,nsteps);
yF=zeros(1,nsteps);
yB=zeros(1,nsteps);
yM=zeros(1,nsteps);
yP=zeros(1,nsteps);
yC=zeros(1,nsteps);
%yy=zeros(1,nsteps);

%Initial Conditions
xi=1;
yF=1;
yB=1;
yM=1;
yP=1;
yC=1;
yy=1

fprintf('\n  Time         Forward_Mehtod          Backward_Method          Modified_Method\')
fprintf('\n  -----------------------------------------------------------------------------\n')

for j=1:nsteps
    xi=j*dx;
    %Forward Method
    yF=yF+(2*(xi)*yF)*dx;
    xinew(j)=xi;
    yFnew(j)=yF;
    %Backward Method
    yB=-yB/(dx*(2*xi)-1);
    xinew(j)=xi;
    yBnew(j)=yB;
    %Modified
    yM=yM+(dx/2)*((2*xi*yM)+2*(xi+dx)*(yM+2*xi*yM*dx));
    xinew(j)=xi;
    yMnew(j)=yM;
    %Predictor Corrector Method
    yP=yP+dx*(2*xi*yP);       %Predictor      
    yC=yC+dx*(2*xi*yP);
    xinew(j)=xi;
    yPnew(j)=yP;
    yCnew(j)=yC;
    %Analytic Solution
    yy=exp((xi^2)-1);

end

fprintf('%6.2g %20.2f %20.2f %20.2f %20.2f %20.2f\n', [xi, yFnew, yBnew, yMnew, yPnew, yCnew, yy])

plot(xinew,yFnew,'b--','linewidth',3)
hold on
plot(xinew,yBnew,'g--','linewidth',3)
plot(xinew,yMnew,'r--','linewidth',3)
plot(xinew,yPnew,'*','linewidth',3)

legend('Forward Method','Backward Method','Modified Method','Predictor-Corrector')
xlabel('time')
ylabel('y')

[Walter Roberson]

What leads you to suspect that there is something that needs to be fixed?

[Romina Doubnia]

The values that I get on the table are not right. I should be getting for Time something like:

Time         Forward_Mehtod          Backward_Method          Modified_Method
-----------------------------------------------------------------------------
 1                 1.00                 1.02                 1.03                 1.05                 1.08
1.1                 1.01
1.15                 1.19                 1.15                 1.31              
1.20                 1.03                 1.05                 1.08                 1.25                 1.15
1.30                1.26                 1.32                 1.01                 1.35                 1.04
1.45                1.09                 1.13                 1.17                 1.50                 1.28

If the time column gets fixed is likely that my values for forward, backward, and modified could show more accurately to the results I did by hand.

[Torsten]

If it helps: The analytical solution for your problem is

y(x) = exp(x^2-1)

Best wishes

Torsten.

[Romina Doubnia]

Thank you very much

Answer 1

[Walter Roberson]

You have

fprintf('%6.2g %20.2f %20.2f %20.2f %20.2f %20.2f\n', [xi, yFnew, yBnew, yMnew, yPnew, yCnew, yy])

xi is a scalar. yFnew and yBnew and yMnew and yPnew and yCnew are vectors of length nsteps. yy is a scalar.

You will want to adjust the xi and yi items to be vectors the same length as the other items. And inside the [] in that call, you will want to use semi-colons instead of commas.

Answer 2

[John D'Errico]

I'm not going to debug that confusing mess of code. Sorry. And it looks as if my statement about the solution being an oscillatory one might make sense, if your predictions were correct. The table that you produce is a complex mess, making it appear as if the solution is oscillatory. In fact, that table is utterly confusing, showing 5 columns, with only 3 labels on top.

In fact, the solution is easy to compute analytically. (Hint: separation of variables.) It is:

   y = exp(x^2-1)

So plotting the analytical solution, we see:

This is what you should hope to see. And, since the problem is not stiff at all, Euler's method should work reasonably well, subject to the caveat that the step size may be insufficiently small.

So, I'll give you a simple code for Euler's method here. (You did make a credible stab at the problem, although I won't write the rest of it for you.)

dx = 0.05;
x = 1:.05:1.5;
y = zeros(size(x));
y(1) = 1;
for n = 2:numel(x);
  y(n) = y(n-1) + 2*x(n-1)*y(n-1)*dx;
end 
ezplot(@(x) exp(x.^2 - 1),[1,1.5])
hold on
plot(x,y,'ro')

That is what you SHOULD get. You can add loops to do the other methods, based on what I did.

I would in fact recommend that you keep the methods separate, since as you wrote it, you are probably screwing things up by trying to do all three methods in one loop.

[Romina Doubnia]

Your criticism was well appreciated it. Thank you. I'll recheck again taking your opinion into consideration.

[Romina Doubnia]

For each loop, I am supposed to have different values, correct? (Just making sure).

[John D'Errico]

Yes. Each method would give you different predictions, because they have subtly different behavior and accuracy. As you can see, the simple Euler method loop I wrote is not that accurate. But the other methods might be more accurate (or potentially less so.) All will probably have the correct general behavior.

Note that if you cut the step size, you would expect the predictions to follow the correct curve more accurately.