Changeset 15451 in project
 Timestamp:
 08/14/09 02:16:43 (10 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

release/4/nemo/trunk/examples/carelli05_run.m
r15166 r15451 1 1 #!/usr/bin/env octave 2 2 3 function ys = euler(f,t,y,h) 4 t0=t(1); tn=t(2); 5 tx=t0; ys=zeros(ceil((tnt0)/h),length(y)); 6 ys(1,:)=y; i=1; 7 while (tx<tn) 8 i = i+1; tx=tx+h; 9 ys(i,:) = ys(i1,:)' + (h * feval(f,tx,ys(i1,:))); 10 endwhile 11 endfunction 3 function 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 % 17 TOL = 0.00001; 18 19 x(1) = t; 20 nstep = ( t(2)  t(1) ) / h; 21 y(:,1) = y_initial; 22 23 for 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 41 end 12 42 13 43 ## Carelli 05 model driver for Octave
Note: See TracChangeset
for help on using the changeset viewer.