/usr/share/octave/site/m/sundialsTB/idas/examples_ser/midasRoberts_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 | function midasRoberts_dns
%midasRoberts_dns - IDAS example problem (serial, dense)
% The following is a simple example problem, with the coding
% needed for its solution by IDAS. 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
% 1 = y1 + y2 + y3
% 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.
% Radu Serban <radu@llnl.gov>
% Copyright (c) 2007, 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];
yp0 = [-0.04;0.04;0.0];
options = IDASetOptions('UserData', data,...
'RelTol',1.e-4,...
'AbsTol',[1.e-8; 1.e-14; 1.e-6],...
'LinearSolver','Dense',...
'JacobianFn',@djacfn);
options = IDASetOptions(options,'RootsFn',@rootfn, 'NumRoots',2);
%mondata.sol = true;
mondata.updt = 100;
options = IDASetOptions(options,'MonitorFn',@IDAMonitor,'MonitorData',mondata);
IDAInit(@resfn,t0,y0,yp0,options);
t1 = 0.4;
tmult = 10.0;
nout = 12;
fprintf('-----------------------------------------------------------------------\n');
fprintf(' t y1 y2 y3');
fprintf(' | nst k h\n');
fprintf('-----------------------------------------------------------------------\n');
iout = 0;
tout = t1;
while iout < nout
[status,t,y] = IDASolve(tout,'Normal');
% Extract statistics
si = IDAGetStats;
% Print output
if(status == 2)
fprintf(' ... Root found %d %d\n',si.RootInfo.roots(1), si.RootInfo.roots(2));
end
fprintf('%10.4e %12.4e %12.4e %12.4e | %3d %1d %12.4e\n',...
t, y(1), y(2), y(3), si.nst, si.qlast, si.hlast);
% Update output time
if(status == 0)
iout = iout+1;
tout = tout*tmult;
end
end
si = IDAGetStats;
IDAFree;
function [rr, flag, new_data] = resfn(t, y, yp, data)
% DAE residual function
r1 = data.p(1);
r2 = data.p(2);
r3 = data.p(3);
rr(1) = -r1*y(1) + r2*y(2)*y(3) - yp(1);
rr(2) = r1*y(1) - r2*y(2)*y(3) - r3*y(2)*y(2) - yp(2);
rr(3) = y(1) + y(2) + y(3) - 1.0;
flag = 0;
new_data = [];
function [J, flag, new_data] = djacfn(t, y, yp, rr, cj, data)
% Dense Jacobian function
r1 = data.p(1);
r2 = data.p(2);
r3 = data.p(3);
J(1,1) = -r1 - cj;
J(2,1) = r1;
J(3,1) = 1.0;
J(1,2) = r2*y(3);
J(2,2) = -r2*y(3) - 2*r3*y(2) - cj;
J(3,2) = 1.0;
J(1,3) = r2*y(2);
J(2,3) = -r2*y(2);
J(3,3) = 1.0;
flag = 0;
new_data = [];
function [g, flag, new_data] = rootfn(t,y,yp,data)
% Root finding function
g(1) = y(1) - 0.0001;
g(2) = y(3) - 0.01;
flag = 0;
new_data = [];
|