1

I am trying to integrate dy3dt = @(t,y) -41*(y - t) + 1; over [0, 1] with h = 0.05 using backward Euler. The initial condition, y(0) is 1. 

yvals3 = [1]; %inital condition
dy3dt = @(t,y) -41*(y - t) + 1;
for t = 0:0.05:0.95; 

    index = round(t*20 + 1); %convert t into index of the last step
    oldy = yvals3(index);
    newy1 = oldy + h.*dy3dt(t, oldy);
    newy2 = oldy + h.*dy3dt(t+ 0.05, newy1); 
    yvals3 = [yvals3; newy2];
end
tvals3 = 0:0.05:1;

When I graph tvals3 and yvals3, I am getting an unusual exponential graph, so my approach is probably incorrect. Any insight on how to implement this using fzero?

Xiao Wang
  • 13
  • 4
  • I don't see the [backward Euler method](http://en.wikipedia.org/wiki/Backward_Euler_method) here. This method involves solving a certain algebraic equation at every step (hence, the mention of `fzero`). Begin by clarifying what that equation is. The code you have is for some form of modified Euler's method. (Second-order Runge-Kutta)/ –  Feb 17 '15 at 19:36
  • I am not sure, that was my initial attempt based on some online sources, which now I realize are incorrect. – Xiao Wang Feb 17 '15 at 19:39
  • Do you see how I can use fzero to get the backward euler? I belief I have to use fzero(newy1 - function, xval). – Xiao Wang Feb 17 '15 at 19:41
  • You need to solve the equation `newy - h*f(x+h,newy) = oldy` for newy. –  Feb 17 '15 at 19:43
  • First newy must be estimated right? I presume we use fzero to that – Xiao Wang Feb 17 '15 at 19:48
  • No, this is not a predictor-corrector method. You use `fzero` to solve for `newy`, and that's it. –  Feb 17 '15 at 19:49
  • This appears to be a problem from UW AMATH 301 HW4... (http://courses.washington.edu/am301/hw/2015win/hw4.pdf). – Anubian Noob Feb 27 '15 at 17:41

1 Answers1

1

The usual (forward) Euler's method can be expressed as going from a known point on a tangent, and getting new point: 

newy = oldy + tstep*dydt(oldt, oldy);

The backward Euler method does everything backwards: it goes from a new (yet unknown) point on a tangent, backward, and hits the old point: 

oldy = newy - tstep*dydt(newt, newy);

This last line is not a command that can be executed, since newy is unknown. Instead, it's an equation to be solved: 

newy = fzero(@(y) y - tstep*dydt(newt,y) - oldy, oldy)

(The second parameter, oldy, is an initial guess to the root, which fzero needs.)

Here's how this could work in your setup. I wrote the code in a way that may not be the most Matlab-efficient, but (I hope) makes the logic of computation clear.

dydt = @(t,y) -41*(y - t) + 1;   % right hand side
tvals = [0];                     % initial t
yvals = [1];                     % initial y
tstep = 0.05;                    % time step 
numsteps = 20;                   % number of steps to go
for i=1:numsteps; 
    oldt = tvals(end);           % pick the last entry of array
    oldy = yvals(end);
    newt = oldt + tstep;         % step forward in t 
    newy = fzero(@(y) y - tstep*dydt(newt,y) - oldy, oldy);  % backward Euler 
    tvals = [tvals; newt];       
    yvals = [yvals; newy];
end
plot(tvals, yvals);