Embedded Files
strangerland®gmail.com
Here we are considering the 1st order linear differential equation
where Φ is the V×E directed incidence matrix of the network (edge orientation is unimportant here), K is an invertible E×E diagonal matrix of edge conductivities, C is an invertible V×V diagonal matrix of vertex capacities, x is a V×1 column vector representing node temperature, and f is a V×1 column vector inhomogeneous term. The diagonal elements of C would be the mass times the heat capacity of each node, and the diagonal components of K would have to be a function of the end point thermal conductivities, the edge cross sectional area, and the geometry of each end point. One possibility is
where vi is the volume of the end node, aj is the cross sectional area of the edge and ki the the thermal conductivity of the end node. Another possibility might be
The differential equation can be rewritten as
or with the change of variables y = Cx
The variable y would represent the "heat" in each node. This expresses a naive relationship between heat and temperature. Such a model could only be applied over a limited temperature range. The heat capacity and thermal conductivity would properly be functions of the node temperature. Note that the matrix
has the requisite properties of a transition rate matrix for a finite state continuous time Markov process.
Output from the following java which uses Runge-Kutta 4th order integration to step forward.
//***************************************************************************************************************************************************************
//ThermalNetwork.java
package com.smkltd.app;
import Jama.Matrix;
public class ThermalNetwork
{
int V;
int E;
Matrix EdgeConductivity;
Matrix VertexCapacity;
Matrix Inhomogeneous;
Matrix CurrentTemp;
Matrix dydt;
TwoBit Incidence;
public ThermalNetwork(int nrows, int ncols)
{
this.V = nrows;
this.E = ncols;
this.EdgeConductivity = new Matrix(new double[1][ncols]);
this.VertexCapacity = new Matrix(new double[nrows][1]);
this.Inhomogeneous = new Matrix(new double[nrows][1]);
this.CurrentTemp = new Matrix(new double[nrows][1]);
this.dydt = new Matrix(new double[nrows][1]);
this.Incidence = new TwoBit(nrows,ncols);
}
public void rk4LDE1(double step)
{
double hh = 0.5*step;
double h6 = step/6;
Matrix dym = new Matrix(new double[V][1]);
Matrix dyt = new Matrix(new double[V][1]);
Matrix yt = new Matrix(new double[V][1]);
dydt = Incidence.LaplacianMod.times(CurrentTemp).plus(Inhomogeneous);
yt = CurrentTemp.plus(dydt.times(hh));
dyt = Incidence.LaplacianMod.times(yt).plus(Inhomogeneous);
yt = CurrentTemp.plus(dyt.times(hh));
dym = Incidence.LaplacianMod.times(yt).plus(Inhomogeneous);
yt = CurrentTemp.plus(dym.times(step));
dym = dym.plus(dyt);
dyt = Incidence.LaplacianMod.times(yt).plus(Inhomogeneous);
CurrentTemp = CurrentTemp.plus(dydt.plus(dyt).plus(dym.times(2.0)).times(h6));
}
public static void main( String[] args )
{
System.out.println(" ");
int[][] MyIncidence = {
{1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
{-1, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
{0, -1, 0, -1, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
{0, 0, -1, 0, -1, 0, -1, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
{0, 0, 0, 0, 0, -1, 0, -1, 0, -1, 0, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0},
{0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, -1, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0},
{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, -1, 0, 0, 0, -1, 0, 1, 0, 0, 0},
{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, -1, 0, 1, 1, 0},
{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, -1, -1, 0, 1},
{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, -1, -1}
};
Matrix MyEC = new Matrix(new double[][] {{1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1}}); // Konductivity of edges
Matrix MyVC = new Matrix(new double[][] {{1,1,1,1,1,1,1,1,1,1}}); // mass times heat capacity of each vertex
Matrix MyInitTemp = new Matrix(new double[][] {{100,150,200,75,60,25,212,125,40,12}});
ThermalNetwork MyThermalNetwork = new ThermalNetwork(MyIncidence.length,MyIncidence[0].length);
for (int row = 0; row < MyThermalNetwork.V; row++)
{
for (int col = 0; col < MyThermalNetwork.E; col++)
{
MyThermalNetwork.Incidence.Set(row,col,MyIncidence[row][col]);
}
}
MyThermalNetwork.EdgeConductivity = MyEC;
MyThermalNetwork.VertexCapacity = MyVC.transpose();
MyThermalNetwork.CurrentTemp = MyInitTemp.transpose().copy();
MyThermalNetwork.Incidence.Laplacian(MyThermalNetwork.EdgeConductivity, MyThermalNetwork.VertexCapacity);
System.out.printf(" ");
MyThermalNetwork.CurrentTemp.transpose().print(8,4);
for (int ii = 0; ii < 125; ii++)
{
MyThermalNetwork.rk4LDE1(.05);
MyThermalNetwork.CurrentTemp.transpose().print(8,4);
}
System.out.println(" ");
System.out.println(" ");
System.out.println(" Printing LaplacianMod = (C^-1)Phi.K.PhiT for MyThermalNetwork");
MyThermalNetwork.Incidence.LaplacianMod.print(6,4);
}
}
//***********************************************************************************************************************
//TwoBit.java
package com.smkltd.app;
import java.util.BitSet; //to make a two bit two's complement array like object for totally unimodular matrix
import Jama.Matrix;
public class TwoBit
{
public int Rows;
public int Cols;
public int[][] Laplacian;
public int[][] Adjacency;
public Matrix LaplacianMod;
public Matrix PKQTC;
private BitSet[][] twobit;
public TwoBit(int nrows, int ncols)
{
this.Rows = nrows;
this.Cols = ncols;
this.twobit = new BitSet[2][nrows];
for (int ii=0;ii<nrows;ii++)
{
this.twobit[0][ii] = new BitSet(ncols);
this.twobit[1][ii] = new BitSet(ncols);
}
}
public int Get(int r, int c)
{
int ii;
if (r>=Rows) {
throw new IllegalArgumentException("row index out of bounds in TwoBit.Get");
} else if (c>=Cols) {
throw new IllegalArgumentException("col index out of bounds in TwoBit.Get");
} else {
if ((twobit[0][r].get(c))&(twobit[1][r].get(c))) {
ii = -1;
} else if (!(twobit[0][r].get(c)) & !(twobit[1][r].get(c))) {
ii=0;
} else if (!(twobit[0][r].get(c)) & (twobit[1][r].get(c))) {
ii=1;
} else {
ii=-2;
}
}
return ii;
}
public void Set(int r, int c, int pm1v0)
{
if ((pm1v0!=-1)&(pm1v0!=0)&(pm1v0!=1)) {
throw new IllegalArgumentException("trying to set value <> -1,0,1 in TwoBit.Set");
} else if (r>=Rows) {
throw new IllegalArgumentException("row index out of bounds in TwoBit.Set");
} else if (c>=Cols) {
throw new IllegalArgumentException("col index out of bounds in TwoBit.Set");
} else {
if (pm1v0 == -1) {
twobit[0][r].set(c);
twobit[1][r].set(c);
} else if (pm1v0 == 0) {
twobit[0][r].clear(c);
twobit[1][r].clear(c);
} else if (pm1v0 == 1) {
twobit[0][r].clear(c);
twobit[1][r].set(c);
}
}
}
public int Pos(int r, int c)
{
int ii;
if (r>=Rows) {
throw new IllegalArgumentException("row index out of bounds in TwoBit.Get");
} else if (c>=Cols) {
throw new IllegalArgumentException("col index out of bounds in TwoBit.Get");
} else {
ii = Get(r,c);
if (ii<0) {
ii=0;
}
}
return ii;
}
public int Neg(int r, int c)
{
int ii;
if (r>=Rows) {
throw new IllegalArgumentException("row index out of bounds in TwoBit.Get");
} else if (c>=Cols) {
throw new IllegalArgumentException("col index out of bounds in TwoBit.Get");
} else {
ii = Get(r,c);
if (ii>0) {
ii=0;
}
}
return -ii;
}
public void Adjacency()
{
Adjacency = new int[Rows][Rows];
for (int ii=0; ii<Rows; ii++)
{
for (int jj=0; jj<Rows; jj++)
{
for (int kk=0; kk<Cols; kk++)
{
Adjacency[ii][jj] = Adjacency[ii][jj] + Pos(ii,kk)*Neg(jj,kk);
}
}
}
}
public void Adjacency(Matrix EdgeWeight, boolean normalized)
{
PKQTC = new Matrix(new double[Rows][Rows]);
for (int ii=0; ii<Rows; ii++)
{
for (int jj=0; jj<Rows; jj++)
{
for (int kk=0; kk<Cols; kk++)
{
PKQTC.set(ii,jj,PKQTC.get(ii,jj) + EdgeWeight.get(0,kk)*Pos(ii,kk)*Neg(jj,kk));
}
}
}
if (normalized) {
double w;
for (int ii=0; ii<Rows; ii++)
{
w=0;
for (int jj=0; jj<Rows; jj++)
{
w = w + PKQTC.get(ii,jj);
}
for (int jj=0; jj<Rows; jj++)
{
PKQTC.set(ii,jj,PKQTC.get(ii,jj)/w);
}
}
}
}
public void Adjacency(Matrix EdgeWeight, Matrix VertexWeight, boolean normalized)
{
Adjacency(EdgeWeight,false);
for (int ii=0; ii<Rows; ii++)
{
PKQTC.set(ii,ii,VertexWeight.get(ii,0));
}
if (normalized) {
double w;
for (int ii=0; ii<Rows; ii++)
{
w=0;
for (int jj=0; jj<Rows; jj++)
{
w = w + PKQTC.get(ii,jj);
}
for (int jj=0; jj<Rows; jj++)
{
PKQTC.set(ii,jj,PKQTC.get(ii,jj)/w);
}
}
}
}
public void Laplacian()
{
Laplacian = new int[Rows][Rows];
for (int ii=0; ii<Rows; ii++)
{
for (int jj=0; jj<Rows; jj++)
{
Laplacian[ii][jj] = Laplacian(ii,jj);
}
}
}
public int Laplacian(int r, int c)
{
//also called the Kirchhoff matrix
int lp;
if (r>=Rows) {
throw new IllegalArgumentException("row index out of bounds in TwoBit.Laplacian");
} else if (c>=Rows) {
throw new IllegalArgumentException("col index out of bounds in TwoBit.Laplacian");
} else {
lp = 0;
for (int ii=0; ii<Cols; ii++)
{
lp = lp + Get(r,ii)*Get(c,ii);
}
}
return lp;
}
public void Laplacian(Matrix K)
{
LaplacianMod = new Matrix(new double[Rows][Rows]);
if ((K.getColumnDimension()!=Cols)|(K.getRowDimension()!=1)) {
throw new IllegalArgumentException("K has incorrect dimensions in TwoBit.Laplacian");
} else {
for (int ii=0; ii<Rows; ii++)
{
for (int jj=0; jj<Rows; jj++)
{
for (int kk=0; kk<Cols; kk++)
{
LaplacianMod.set(ii,jj,LaplacianMod.get(ii,jj)-Get(ii,kk)*Get(jj,kk)*K.get(0,kk));
}
}
}
}
}
public void Laplacian(Matrix K, Matrix C)
{
if ((C.getColumnDimension()!=1)|(C.getRowDimension()!=Rows)) {
throw new IllegalArgumentException("C has incorrect dimensions in TwoBit.Laplacian");
} else {
Laplacian(K);
for (int ii=0; ii<Rows; ii++)
{
for (int jj=0; jj<Rows; jj++)
{
LaplacianMod.set(ii,jj,LaplacianMod.get(ii,jj)/C.get(ii,0));
}
}
}
}
public void xchange(int pvtr, int pvtc)
{
if (pvtr>=Rows) {
throw new IllegalArgumentException("row index out of bounds in TwoBit.xchange()");
} else if (pvtc>=Cols) {
throw new IllegalArgumentException("col index out of bounds in TwoBit.xchange()");
} else if (Get(pvtr,pvtc)==0) {
throw new IllegalArgumentException("Illegal pivot element in TwoBit.xchange()");
} else {
for (int ii=0; ii<Rows; ii++) //need to catch a not totally unimodular condition though set() will bomb in this case
{
for (int jj=0; jj<Cols; jj++)
{
if ((ii!=pvtr) & (jj!=pvtc)) {
Set(ii,jj,(Get(ii,jj)-Get(ii,pvtc)*Get(pvtr,jj)/Get(pvtr,pvtc)));
}
}
if (ii!=pvtr) {
Set(ii,pvtc,(Get(ii,pvtc)/Get(pvtr,pvtc)));
}
}
for(int jj=0; jj<Cols; jj++)
{
if (jj!=pvtc) {
Set(pvtr,jj,(-Get(pvtr,jj)/Get(pvtr,pvtc)));
}
}
Set(pvtr,pvtc,(1/Get(pvtr,pvtc)));
}
}
public static void main( String[] args )
{
int v = 5;
int e;
e=(int)(v*(v-1)/2);
TwoBit mytwobit = new TwoBit(v,e);
//mytwobit.Create(v,e);
int kk=-1;
for (int ii=0; ii<v-1;ii++)
{
for (int jj=ii+1; jj<v;jj++)
{
kk=kk+1;
mytwobit.Set(ii,kk,1);
mytwobit.Set(jj,kk,-1);
}
}
System.out.println(" ");
for (int ii=0; ii<mytwobit.Rows;ii++)
{
for (int jj=0; jj<mytwobit.Cols;jj++)
{
System.out.printf("%2d ", mytwobit.Get(ii,jj));
}
System.out.println("");
}
}
}
Page updated
Google Sites
Report abuse