/usr/share/octave/site/m/sundialsTB/cvodes/examples_ser/mcvsRoberts_dns.m is in octave-sundials 2.5.0-3+b1.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 | function mcvsRoberts_dns()
%mcvsRoberts_dns - CVODES example problem (serial, dense)
% The following is a simple example problem, with the coding
% needed for its solution by CVODES. The problem is from
% chemical kinetics, and consists of the following three rate
% equations:
% dy1/dt = -.04*y1 + 1.e4*y2*y3
% dy2/dt = .04*y1 - 1.e4*y2*y3 - 3.e7*(y2)^2
% dy3/dt = 3.e7*(y2)^2
% on the interval from t = 0.0 to t = 4.e10, with initial
% conditions: y1 = 1.0, y2 = y3 = 0. The problem is stiff.
% While integrating the system, we also use the rootfinding
% feature to find the points at which y1 = 1e-4 or at which
% y3 = 0.01. This program solves the problem with the BDF method,
% Newton iteration with the CVDENSE dense linear solver, and a
% user-supplied Jacobian routine. It uses a scalar relative
% tolerance and a vector absolute tolerance and also illustrates
% the rootfinding capability in CVODES.
%
% Output is printed in decades from t = .4 to t = 4.e10.
% Run statistics (optional outputs) are printed at the end.
% Radu Serban <radu@llnl.gov>
% Copyright (c) 2005, The Regents of the University of California.
% $Revision: 1.1 $Date: 2007/10/26 16:30:48 $
data.p = [0.04; 1.0e4; 3.0e7];
t0 = 0.0;
y0 = [1.0;0.0;0.0];
options = CVodeSetOptions('UserData', data,...
'RelTol',1.e-8,...
'AbsTol',[1.e-8; 1.e-14; 1.e-6],...
'LinearSolver','Dense',...
'JacobianFn',@djacfn,...
'RootsFn',@rootfn, 'NumRoots',2);
mondata.sol = true;
mondata.mode = 'text';
mondata.skip = 9;
mondata.updt = 100;
options = CVodeSetOptions(options,'MonitorFn',@CVodeMonitor,'MonitorData',mondata);
CVodeInit(@rhsfn, 'BDF', 'Newton', t0, y0, options);
t1 = 0.4;
tmult = 10.0;
nout = 12;
iout = 0;
tout = t1;
while iout < nout
[status,t,y] = CVode(tout,'Normal');
% Extract statistics
si = CVodeGetStats;
% Print output
fprintf('t = %0.2e order = %1d step = %0.2e',t, si.qlast, si.hlast);
if(status == 2)
fprintf(' ... Root found %d %d\n',si.RootInfo.roots(1), si.RootInfo.roots(2));
else
fprintf('\n');
end
fprintf('solution = [ %14.6e %14.6e %14.6e ]\n\n', y(1), y(2), y(3));
% Update output time
if(status == 0)
iout = iout+1;
tout = tout*tmult;
end
end
si = CVodeGetStats;
CVodeFree;
% ===========================================================================
function [yd, flag, new_data] = rhsfn(t, y, data)
% Right-hand side function
r1 = data.p(1);
r2 = data.p(2);
r3 = data.p(3);
yd(1) = -r1*y(1) + r2*y(2)*y(3);
yd(3) = r3*y(2)*y(2);
yd(2) = -yd(1) - yd(3);
flag = 0;
new_data = [];
return
% ===========================================================================
function [J, flag, new_data] = djacfn(t, y, fy, data)
% Dense Jacobian function
r1 = data.p(1);
r2 = data.p(2);
r3 = data.p(3);
J(1,1) = -r1;
J(1,2) = r2*y(3);
J(1,3) = r2*y(2);
J(2,1) = r1;
J(2,2) = -r2*y(3) - 2*r3*y(2);
J(2,3) = -r2*y(2);
J(3,2) = 2*r3*y(2);
flag = 0;
new_data = [];
return
% ===========================================================================
function [g, flag, new_data] = rootfn(t,y,data)
% Root finding function
g(1) = y(1) - 0.0001;
g(2) = y(3) - 0.01;
flag = 0;
new_data = [];
|