#Maple worksheet with Lectures, Part 17 # # follows section 7.10 of Yeargers # except uses symbolic definition for # the matrix A and then substitutes # constants # ### first just show equations 7.10.1diff(x1(t),t)=-(a01+a21+a31)*x1(t)+a12*x2(t)+a13*x3(t)+Il(t);

NiMvLUklZGlmZkdJKnByb3RlY3RlZEdGJjYkLUkjeDFHNiI2I0kidEdGKkYsLCoqJiwoSSRhMDFHRioiIiJJJGEyMUdGKkYxSSRhMzFHRipGMUYxRihGMSEiIiomSSRhMTJHRipGMS1JI3gyR0YqRitGMUYxKiZJJGExM0dGKkYxLUkjeDNHRipGK0YxRjEtSSNJbEdGKkYrRjE=

diff(x2(t),t)=a21*x1(t)-(a02+a12+a32)*x2(t)+a23*x3(t);

NiMvLUklZGlmZkdJKnByb3RlY3RlZEdGJjYkLUkjeDJHNiI2I0kidEdGKkYsLCgqJkkkYTIxR0YqIiIiLUkjeDFHRipGK0YwRjAqJiwoSSRhMDJHRipGMEkkYTEyR0YqRjBJJGEzMkdGKkYwRjBGKEYwISIiKiZJJGEyM0dGKkYwLUkjeDNHRipGK0YwRjA=

diff(x3(t),t)=a31*x1(t)+a32*x2(t)-(a03+a13+a23)*x3(t);

NiMvLUklZGlmZkdJKnByb3RlY3RlZEdGJjYkLUkjeDNHNiI2I0kidEdGKkYsLCgqJkkkYTMxR0YqIiIiLUkjeDFHRipGK0YwRjAqJkkkYTMyR0YqRjAtSSN4MkdGKkYrRjBGMComLChJJGEwM0dGKkYwSSRhMTNHRipGMEkkYTIzR0YqRjBGMEYoRjAhIiI=

# define a generic three-compartment # first order model Asymb:=matrix(3,3,[-(a01+a21+a31),a12,a13,a21,-(a02+a12+a32),a23,a31,a32,-(a03+a13+a23)]);

NiM+SSZBc3ltYkc2Ii1JJ21hdHJpeEc2JEkqcHJvdGVjdGVkR0YpSShfc3lzbGliR0YlNiM3JTclLChJJGEwMUdGJSEiIkkkYTIxR0YlRjBJJGEzMUdGJUYwSSRhMTJHRiVJJGExM0dGJTclRjEsKEkkYTAyR0YlRjBGM0YwSSRhMzJHRiVGMEkkYTIzR0YlNyVGMkY4LChJJGEwM0dGJUYwRjRGMEY5RjA=

Xsymb:=[x1(t),x2(t),x3(t)];

NiM+SSZYc3ltYkc2IjclLUkjeDFHRiU2I0kidEdGJS1JI3gyR0YlRiktSSN4M0dGJUYp

evalm(Asymb&*Xsymb);

NiMtSSd2ZWN0b3JHNiI2IzclLCgqJiwoSSRhMDFHRiUhIiJJJGEyMUdGJUYsSSRhMzFHRiVGLCIiIi1JI3gxR0YlNiNJInRHRiVGL0YvKiZJJGExMkdGJUYvLUkjeDJHRiVGMkYvRi8qJkkkYTEzR0YlRi8tSSN4M0dGJUYyRi9GLywoKiZGLUYvRjBGL0YvKiYsKEkkYTAyR0YlRixGNUYsSSRhMzJHRiVGLEYvRjZGL0YvKiZJJGEyM0dGJUYvRjpGL0YvLCgqJkYuRi9GMEYvRi8qJkZBRi9GNkYvRi8qJiwoSSRhMDNHRiVGLEY5RixGQ0YsRi9GOkYvRi8=

# define the rates as given in # Table 7.10.1 rates:={a01=0.0211,a02=0.0162,a03=0,a12=0.0124,a13=0.000035,a21=0.0111,a23=0,a31=0.0039,a32=0}

NiM+SSZyYXRlc0c2IjwrL0kkYTAxR0YlJCIkNiMhIiUvSSRhMDJHRiUkIiRpIkYrL0kkYTAzR0YlIiIhL0kkYTEyR0YlJCIkQyJGKy9JJGExM0dGJSQiI04hIicvSSRhMjFHRiUkIiQ2IkYrL0kkYTIzR0YlRjIvSSRhMzFHRiUkIiNSRisvSSRhMzJHRiVGMg==

;# substitute the numeric values # "subs" does not work on matrices # so we need to use map2 to perform # the substitution on each element A:=map2(subs,rates,Asymb);

NiM+SSJBRzYiLUknbWF0cml4RzYkSSpwcm90ZWN0ZWRHRilJKF9zeXNsaWJHRiU2IzclNyUkISRoJCEiJSQiJEMiRjAkIiNOISInNyUkIiQ2IkYwJCEkJ0dGMCIiITclJCIjUkYwRjskISNORjU=

with(linalg):

Warning, the protected names norm and trace have been redefined and unprotected

# from here follow syntax on page 223 etA:=exponential(A,t);

NiM+SSRldEFHNiItSSdtYXRyaXhHNiRJKnByb3RlY3RlZEdGKUkoX3N5c2xpYkdGJTYjNyU3JSwoLUkkZXhwR0YoNiMsJEkidEdGJSQhK2peKG9ZJSEjNiQiK3Z6ZUFsISM1LUYwNiMsJEYzJCErOTtjLj9GNiQiK1BZLHdNRjktRjA2IywkRjMkIStZRUFqSSEjOSQiK1kkKlEoUiIhIzgsKEZBJCInM2xnRjlGOiQiK2szeEtdRjlGLyQhK3R0UExdRjksKEZBJCIrVSFmKD42ISM3RjokIStHI1FHMydGSUYvJCErJT5fWjYmRkk3JSwoRjokIitnOzkwWEY5Ri8kISslZSVvMFhGOUZBJCInREhhRjksKEYvJCIrNV4ncFokRjlGQSQiJ1pjQkY5RjokIitUIyp6QWxGOSwoRkEkIit4d2RdVkZJRjokIStJXXYkKXlGSUYvJCIrYHQ8TE5GSTclLChGOiQhKnU+IXluRjlGLyQhKkkmSCpwJkY5RkEkIisvOnRaN0Y5LChGQSQiK1JgYTphRjZGOiQhK2QzZzgpKkY2Ri8kIis5YjApUiVGNiwoRi8kIiYjcFchIipGQSQiK3FwTykqKipGOUY6JCIrW0w2Jz0iRkk=

AinvF:=evalm(inverse(A)&*vector([49.3,0,0]));

NiM+SSZBaW52Rkc2Ii1JJ3ZlY3Rvckc2JEkqcHJvdGVjdGVkR0YpSShfc3lzbGliR0YlNiM3JSQhK0dxNCs9ISInJCEraTBSJylwISIoJCErLkMjZSsjISIl

u:=evalm(-AinvF+etA&*AinvF);

NiM+SSJ1RzYiLUkndmVjdG9yRzYkSSpwcm90ZWN0ZWRHRilJKF9zeXNsaWJHRiU2IzclLCokIitHcTQrPSEiJyIiIi1JJGV4cEdGKDYjLCRJInRHRiUkIStqXihvWSUhIzYkIStqKFspKT4oISIoLUYzNiMsJEY2JCErOTtjLj9GOSQhKyVlV0piKUY8LUYzNiMsJEY2JCErWUVBakkhIzkkIStRcCgqW0FGPCwqJCIraTBSJylwRjxGMUY9JCErPExhMzZGMEYyJCIreEokRyhcRjxGRCQhK3FjIXp0KSEiKSwqJCIrLkMjZSsjISIlRjFGPSQiK0lEInltIkY8RjIkIitIISo9IUgnRlVGRCQhK04jPiIzP0ZZ

plot({u[1](t),u[2](t),u[3](t)},t=0..365,thickness=3);

LSUlUExPVEc2Jy0lJ0NVUlZFU0c2JjdaNyQkIiIhRiskISIjISIoNyQkIis9LioqKSk+ISIqJCIqSTNHWSpGLjckJCIrTzEpeihSRjIkIitGJ3p1Iz1GLjckJCIrYjQocCdmRjIkIitQJUgmW0VGLjckJCIrcjcnZiZ6RjIkIitxdCdSVCRGLjckJCIrJjQ9PDkiISIpJCIrISllZ0RZRi43JCQiK2orJXlbIkZHJCIrVyFwNHEmRi43JCQiKzsvNHg9RkckIitRdGhvbkYuNyQkIitwMk1tQUZHJCIrcWQ+MXhGLjckJCIrRTg8ZUVGRyQiK186ZE8mKUYuNyQkIiskKT0rXUlGRyQiK1ZfYW8jKkYuNyQkIitRLCgqUk1GRyQiK3FBUjcqKkYuNyQkIisjUlEqSFFGRyQiKzUzSFs1ISInNyQkIit2JFJJYiVGRyQiK0JiIXk4IkZnbzckJCIrZiZvPEkmRkckIitBXFs3N0ZnbzckJCIrYGA1d2dGRyQiK24pUlpGIkZnbzckJCIrUSlleiVvRkckIisnPnpcSyJGZ283JCQiK0BdIT5rKEZHJCIrJFxAcU8iRmdvNyQkIiteS0BUJClGRyQiKyhIOndSIkZnbzckJCIrKG9vJUciKkZHJCIrV3VGRTlGZ283JCQiK3BtJio9KipGRyQiKyFceitYIkZnbzckJCIrSlEybzVGLiQiKzQiKT5wOUZnbzckJCIrLjREUDZGLiQiK0A8KlFbIkZnbzckJCIrSyI0Jj43Ri4kIitYc2opXCJGZ283JCQiK1RFPipHIkYuJCIrJ0gnPjQ6RmdvNyQkIitINkRxOEYuJCIrJFJHJz46RmdvNyQkIitEOytVOUYuJCIrdG1ZRjpGZ283JCQiKyJIQjJfIkYuJCIrdUghW2AiRmdvNyQkIitNXG8mZiJGLiQiKy4zdlM6RmdvNyQkIitWJioqUW4iRi4kIitobS9ZOkZnbzckJCIrclxzWDxGLiQiK0pZQF06RmdvNyQkIispUik+Qj1GLiQiK0FqMmE6RmdvNyQkIis6Om4uPkYuJCIrLGBcZDpGZ283JCQiKytPc3Q+Ri4kIitXSDBnOkZnbzckJCIrcjxRXD9GLiQiK1s5V2k6RmdvNyQkIitsU2FGQEYuJCIrMShmWGMiRmdvNyQkIitAMSwvQUYuJCIrKXlRamMiRmdvNyQkIitYZip6RiNGLiQiKyk0Q3ljIkZnbzckJCIrTVM5Z0JGLiQiK0g2Q3A6RmdvNyQkIitYeCZSViNGLiQiK11oTHE6RmdvNyQkIitQJ3BGXiNGLiQiK21rTXI6RmdvNyQkIitEXD0lZSNGLiQiKyo9U0BkIkZnbzckJCIrKXBnQW0jRi4kIitoVipHZCJGZ283JCQiKzpSc05GRi4kIitTMl50OkZnbzckJCIrWFReN0dGLiQiK2FAMnU6RmdvNyQkIitKMWYoKUdGLiQiK0wxYnU6RmdvNyQkIithKCk9bUhGLiQiK2B3KVxkIkZnbzckJCIrKTMpKT0vJEYuJCIra2JOdjpGZ283JCQiKyFmLCQ+SkYuJCIrbGRvdjpGZ283JCQiK3FTMic+JEYuJCIrMVEoZmQiRmdvNyQkIiswK2ltS0YuJCIrNSgzaWQiRmdvNyQkIis0T1paTEYuJCIrRHlXdzpGZ283JCQiK04rej5NRi4kIitJI1FtZCJGZ283JCQiK0hXKm9cJEYuJCIrPy4jb2QiRmdvNyQkIitzSnBxTkYuJCIrIz14cGQiRmdvNyQkIiRsJEYrJCIraSVIcmQiRmdvN1k3JEYqJCIjP0ZHNyQkIishZl5cJSoqISM1JCIpa3pbRUZHNyRGMCQiKiFRO1A1Rkc3JCQiK3hhWyQpSEYyJCIqPHJYRyNGRzckRjYkIip5Q2soUkZHNyRGOyQiKi1UJHkmKUZHNyRGQCQiKyZwJ29pOUZHNyRGRSQiK2xoSy5HRkc3JEZLJCIrKSpvI1xWJUZHNyRGVSQiKzJZRDApKUZHNyRGaW4kIitXRCYqcDhGLjckRmNvJCIrIltBZyc9Ri43JEZpbyQiK0c5QzdCRi43JEZecCQiKyoqPnBaRkYuNyRGY3AkIitmdVRpSkYuNyRGaHAkIiswPG9PTkYuNyRGXXEkIitOVGkhKVFGLjckRmJxJCIrKikzT11URi43JEZncSQiKy83Pj5XRi43JEZcciQiKzE5MmJZRi43JEZhciQiK2J3RWBbRi43JEZmciQiKyQzajUsJkYuNyRGW3MkIitoYUl1XkYuNyRGYHMkIitoUkMlSCZGLjckRmVzJCIrdWxCOmFGLjckRmpzJCIrd0h3MmJGLjckRl90JCIrQ2BmJmYmRi43JEZkdCQiK2A3bG5jRi43JEZpdCQiKzJ4VEtkRi43JEZedSQiK0oicFB5JkYuNyRGY3UkIis5Jj47JGVGLjckRmh1JCIrbnA8dWVGLjckRl12JCIrczU2MWZGLjckRmJ2JCIrJkgqKWYkZkYuNyRGZ3YkIitWPl5pZkYuNyRGXHckIisoUSd5JSlmRi43JEZhdyQiKyE+a0wrJ0YuNyRGZnckIitGRDBAZ0YuNyRGW3gkIit3MG9NZ0YuNyRGYHgkIitpYD9aZ0YuNyRGZXgkIitkdCpwMCdGLjckRmp4JCIrTW9DbWdGLjckRl95JCIrL0V2dGdGLjckRmR5JCIrbkJgITMnRi43JEZpeSQiKzRVRCczJ0YuNyRGXnokIitnMVUiNCdGLjckRmN6JCIrZkdyJjQnRi43JEZoeiQiK2R4XSo0J0YuNyRGXVtsJCIrKSo9dy1oRi43JEZiW2wkIitfc08waEYuNyRGZ1tsJCIrS1InejUnRi43JEZcXGwkIitpSCkqNGhGLjckRmFcbCQiKzVdJz02J0YuNyRGZlxsJCIrXTRXOGhGLjckRltdbCQiK0FQI1w2J0YuN1M3JEYqRio3JEZAJCImYGEmISIlNyRGSyQiJyNweiJGX2hsNyRGVSQiJ2xXUUZfaGw3JEZpbiQiJ0Vaa0ZfaGw3JEZjbyQiJ1xjJSpGX2hsNyRGaW8kIih1U0QiRl9obDckRl5wJCIoVHNmIkZfaGw3JEZjcCQiKGJFKD5GX2hsNyRGaHAkIig9T08jRl9obDckRl1xJCIoLCp6RkZfaGw3JEZicSQiKChIY0pGX2hsNyRGZ3EkIig8ISplJEZfaGw3JEZcciQiKGE5LiVGX2hsNyRGYXIkIihkU1klRl9obDckRmZyJCIoUTgnW0ZfaGw3JEZbcyQiKF0kUWBGX2hsNyRGYHMkIihbZHUmRl9obDckRmVzJCIoWUdBJ0ZfaGw3JEZqcyQiKCllWm1GX2hsNyRGX3QkIigvZTYoRl9obDckRmR0JCIocU1jKEZfaGw3JEZpdCQiKF1ALilGX2hsNyRGXnUkIihfUFkpRl9obDckRmN1JCIoJlFJKilGX2hsNyRGaHUkIigjNDslKkZfaGw3JEZddiQiKCwnUikqRl9obDckRmJ2JCIpPHdINUZfaGw3JEZndiQiKV04eDVGX2hsNyRGXHckIillX0I2Rl9obDckRmF3JCIpbVdvNkZfaGw3JEZmdyQiKWtOPTdGX2hsNyRGW3gkIilqQWo3Rl9obDckRmB4JCIpUjo2OEZfaGw3JEZleCQiKXJmYThGX2hsNyRGangkIilLNS05Rl9obDckRl95JCIpLyJvVyJGX2hsNyRGZHkkIil4YSRcIkZfaGw3JEZpeSQiKV1DUjpGX2hsNyRGXnokIilwMyhlIkZfaGw3JEZjeiQiKVI7TDtGX2hsNyRGaHokIilARyFvIkZfaGw3JEZdW2wkIiltK0Y8Rl9obDckRmJbbCQiKXkkKnA8Rl9obDckRmdbbCQiKWs4Pj1GX2hsNyRGXFxsJCIpVjhqPUZfaGw3JEZhXGwkIikmUSsiPkZfaGw3JEZmXGwkIildI1wmPkZfaGw3JEZbXWwkIilNOi4/Rl9obC0lJkNPTE9SRzYsJSRSR0JHJCIjNSEiIiRGK0ZjYW1GZGFtRmRhbUZhYW1GZGFtRmFhbUZhYW1GZGFtLSUrQVhFU0xBQkVMU0c2JFEidDYiUSFGaWFtLSUlRk9OVEc2JCUqSEVMVkVUSUNBR0ZiYW0tJSVWSUVXRzYkO0ZkYW1GW11sOyQhMisrKy8jb0kxUyEjOiQiLi8hb2tAVj9GMi0lKlRISUNLTkVTU0c2IyIiJA==