/usr/share/octave/site/m/sundialsTB/idas/examples_ser/midasReInit_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 | function [] = midasReInit_dns()
%midasReInit_dns - Illustration of the IDAS reinitialization function
% 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 $
fprintf('Example for integrating over a discontinuity in states\n');
fprintf('using the IDAS re-initialization function\n\n');
fprintf('Integrate over t = [ 0 1.0 ] the DAE:');
fprintf(' y1'' + y1 - y2 = 0\n');
fprintf(' y1 + y2 = 0\n');
fprintf('with initial conditions:\n');
fprintf(' y1(0) = 1.0\n');
fprintf(' y2(0) = -1.0\n');
fprintf('until y2(t*) = -0.5. At t*, perturb:\n');
fprintf(' y1(t*) <- y1(t*) - 0.25\n');
fprintf(' y2(t*) <- y2(t*) + 0.25\n');
fprintf('and continue the integration to t = 1.0\n\n');
t0 = 0.0;
tout = 1.0;
y0 = [1.0;-1.0];
yp0 = [-2.0;0.0];
% Set optional inputs
options = IDASetOptions('RelTol',1.e-4,...
'AbsTol',1.e-5,...
'LinearSolver','Dense');
options = IDASetOptions(options,'RootsFn',@my_rootfct, 'NumRoots',1);
% Initialize solver
IDAInit(@my_resfct,t0,y0,yp0,options);
% Initialize arrays
tt = [];
yy1 = [];
yy2 = [];
% Integrate DAE until root is found
t = t0;
while t<tout
[status, t, y] = IDASolve(tout,'OneStep');
tt = [tt;t];
yy1 = [yy1;y(1)];
yy2 = [yy2;y(2)];
if status == 2
break;
end
end
fprintf('');
% Get yp at current time
yp = IDAGet('DerivSolution',t,0);
% Add discontinuity in solution
% (must be consistent with the algebraic constraint)
t0 = t;
y0 = [y(1)-0.25; y(2)+0.25];
yp0 = [yp(1)+0.25; 0.0];
% Reinitialize solver
IDAReInit(t0,y0,yp0,options);
% Integrate to final time
t = t0;
while t<tout
[status, t, y] = IDASolve(tout,'OneStep');
tt = [tt;t];
yy1 = [yy1;y(1)];
yy2 = [yy2;y(2)];
end
% Free memory
IDAFree;
% Plot solution
hp1 = plot(tt,yy1);
hold on
hp2 = plot(tt,yy2,'r');
set(gca,'XLim',[tt(1) tt(end)]);
plot([t0 t0],get(gca,'YLim'),'k:');
legend([hp1,hp2],'y_1','y_2');
%
% ========================================================
%
function [rr, flag] = my_resfct(t, y, yp)
rr(1) = yp(1) + y(1) - y(2);
rr(2) = y(1) + y(2);
flag = 0;
%
% ========================================================
%
function [g, flag] = my_rootfct(t, y, yp)
g(1) = y(2) + 0.5;
flag = 0;
|