Changeset 15451 in project


Ignore:
Timestamp:
08/14/09 02:16:43 (10 years ago)
Author:
Ivan Raikov
Message:

added backward euler method to nemo examples

File:
1 edited

Legend:

Unmodified
Added
Removed
  • release/4/nemo/trunk/examples/carelli05_run.m

    r15166 r15451  
    11#!/usr/bin/env octave
    22
    3 function ys = euler(f,t,y,h)
    4   t0=t(1); tn=t(2);
    5   tx=t0; ys=zeros(ceil((tn-t0)/h),length(y));
    6   ys(1,:)=y; i=1;
    7   while (tx<tn)
    8     i = i+1; tx=tx+h;
    9     ys(i,:) = ys(i-1,:)' + (h * feval(f,tx,ys(i-1,:)));
    10   endwhile
    11 endfunction
     3function ys = be_newton ( f, fy, t, y_initial, h )
     4%
     5%  function [ x, y ] = be_newton ( f, fy, t, y_initial, h )
     6%
     7%  BE_NEWTON uses NSTEP steps of the backward Euler Newton method
     8%  to estimate Y, the solution of an ODE, at the equally spaced points X in
     9%  the range X_RANGE(1) to X_RANGE(2). 
     10%
     11%  Euler's forward method is used as a predictor,
     12%  and Newton's method is used to seek a solution of Euler's backward formula.
     13%
     14%  The name of the derivative function is F.
     15%  The name of the partial derivative function is FY.
     16%
     17TOL = 0.00001;
     18
     19x(1) = t;
     20nstep = ( t(2) - t(1) ) / h;
     21y(:,1) = y_initial;
     22
     23for i = 1 : nstep
     24
     25  x(i+1) = x(i) + h;
     26  yp = y(:,i) + h * feval ( f, x(i), y(:,i) );
     27
     28  it = 0;
     29  resmax = TOL + 1.0;
     30
     31  while ( resmax > TOL & it < 10 )
     32    r = yp - ( y(:,i) + h * feval ( f, x(i+1), yp ) );
     33    rprime = 1.0 - h * feval ( fy, x(i+1), yp );
     34    yp = yp - r / rprime;
     35    resmax = max ( r );
     36    it = it + 1;
     37  end
     38
     39  ys(:,i+1) = yp;
     40
     41end
    1242
    1343## Carelli 05 model driver for Octave
Note: See TracChangeset for help on using the changeset viewer.