PLEIADES_J - Jacobian function for the PLEIADES example problem.

   see also: pleiades, CVDenseJacFn

0001 function [J, flag, new_data] = pleiades_J(t, y, fy, data)
0002 %PLEIADES_J - Jacobian function for the PLEIADES example problem.
0003 %
0004 %   see also: pleiades, CVDenseJacFn
0005 
0006 % Radu Serban <radu@llnl.gov>
0007 % Copyright (c) 2005, The Regents of the University of California.
0008 % $Revision: 1.2 $Date: 2006/02/13 23:01:27 $
0009 
0010 neq = 28;
0011 
0012 J = zeros(neq,neq);
0013 for i = 1:14
0014   J(i,14+i)=1.0;
0015 end
0016 for i = 2:7
0017   mi=i;
0018   for j = 1:i-1
0019     mj = j;
0020     rij = (y(i)-y(j))^2+(y(i+7)-y(j+7))^2;
0021     rij32 = rij^(3/2);
0022     rij52 = rij^(5/2);
0023     fjh = (1.0-3.0*(y(j)-y(i))^2/rij)/rij32;
0024     J(i+14,j)   = mj*fjh;
0025     J(j+14,i)   = mi*fjh;
0026     fjh = (1.0-3.0*(y(j+7)-y(i+7))^2/rij)/rij32;
0027     J(i+21,j+7) = mj*fjh;
0028     J(j+21,i+7) = mi*fjh;
0029     fjh = -3.0*(y(j)-y(i))*(y(j+7)-y(i+7))/rij52;
0030     J(i+14,j+7) = mj*fjh;
0031     J(j+14,i+7) = mi*fjh;
0032     J(i+21,j)   = mj*fjh;
0033     J(j+21,i)   = mi*fjh;
0034   end
0035 end
0036 for i = 1:7
0037   sumxx = 0.0;
0038   sumxy = 0.0;
0039   sumyy = 0.0;
0040   for j = 1:7
0041     if j ~= i
0042       sumxx = sumxx + J(i+14,j);
0043       sumxy = sumxy + J(i+14,j+7);
0044       sumyy = sumyy + J(i+21,j+7);
0045     end
0046   end
0047   J(i+14,i)   = -sumxx;
0048   J(i+14,i+7) = -sumxy;
0049   J(i+21,i)   = -sumxy;
0050   J(i+21,i+7) = -sumyy;
0051 end
0052 
0053 flag = 0;
0054 new_data = [];