Systems of Differential Equations

/*--------------------------------------------------------------

                    SAS Sample Library

        Name: modex10.sas
 Description: Example program from SAS/ETS User's Guide,
              The MODEL Procedure
       Title: Systems of Differential Equations
     Product: SAS/ETS Software
        Keys: nonlinear simultaneous equation models
        PROC: MODEL
       Notes:

--------------------------------------------------------------*/

data fit;
   length type $ 8;
   input type $ time f ;
datalines;
assoc      0         95.3156
assoc      10        97.8261
assoc      20        102.851
assoc      30        101.506
assoc      40        102.737
assoc      50        103.932
assoc      60        105.481
assoc      70        106.458
assoc      80        107.287
assoc      90        106.778
assoc      100       108.601
assoc      110       111.473
assoc      120       115.147
assoc      130       111.667
assoc      140       114.835
assoc      150       113.911
assoc      160       113.759
assoc      170       116.336
assoc      180       116.574
assoc      190       116.567
assoc      200       117.409
assoc      210       121.065
assoc      220       119.332
assoc      230       120.39
assoc      240       126.192
assoc      250       122.187
assoc      260       122.927
assoc      270       123.114
assoc      280       124.94
assoc      290       127.959
assoc      300       129.619
assoc      310       126.638
assoc      320       126.27
assoc      330       128.662
assoc      340       129.566
assoc      350       131.738
assoc      360       131.98
assoc      370       131.361
assoc      380       134.137
assoc      390       133.676
assoc      400       132.277
assoc      410       133.478
assoc      420       135.748
assoc      430       135.881
assoc      440       134.653
assoc      450       135.946
assoc      460       137.99
assoc      470       137.034
assoc      480       140.372
assoc      490       139.477
assoc      500       139.563
assoc      510       140.564
assoc      520       142.086
assoc      530       140.779
assoc      540       140.389
assoc      550       144.758
assoc      560       148.054
assoc      570       142.608
assoc      580       146.63
assoc      590       147.299
assoc      600       144.674
assoc      610       147.017
assoc      620       148.054
assoc      630       148.512
assoc      640       147.753
assoc      650       151.589
assoc      660       150.636
assoc      670       151.332
assoc      680       150.099
assoc      690       150.497
assoc      700       150.895
assoc      710       151.705
assoc      720       153.556
assoc      730       151.709
assoc      740       157.416
assoc      750       153.535
assoc      760       152.944
assoc      770       153.816
assoc      780       154.755
assoc      790       158.671
assoc      800       155.613
assoc      810       155.237
assoc      820       157.842
assoc      830       156.842
assoc      840       157.126
assoc      850       158.536
assoc      860       158.95
assoc      870       159.514
assoc      880       158.838
assoc      890       162.257
assoc      900       162.525
assoc      910       159.284
assoc      920       164.124
assoc      930       161.183
assoc      940       162.96
assoc      950       162.346
assoc      960       171.145
assoc      970       163.673
assoc      980       165.25
assoc      990       164.323
assoc      1000      164.919
assoc      1010      166.521
assoc      1020      164.93
assoc      1030      166.189
assoc      1040      174.801
assoc      1050      168.71
assoc      1060      167.583
assoc      1070      172.601
assoc      1080      166.945
assoc      1090      168.205
assoc      1100      170.792
assoc      1110      174.314
assoc      1120      169.344
assoc      1130      170.462
assoc      1140      171.963
assoc      1150      172.162
assoc      1160      170.211
assoc      1170      171.167
assoc      1180      170.612
assoc      1190      172.059
assoc      1200      176.471
assoc      1210      173.598
assoc      1220      173.896
assoc      1230      175.865
assoc      1240      176.173
assoc      1250      175.497
assoc      1260      176.433
assoc      1270      172.532
assoc      1280      175.219
assoc      1290      175.07
assoc      1300      174.613
assoc      1310      174.751
assoc      1320      175.179
assoc      1330      176.77
assoc      1340      181.655
assoc      1350      178.105
assoc      1360      174.831
assoc      1370      176.513
assoc      1380      178.623
assoc      1390      176.737
assoc      1400      177.368
assoc      1410      178.847
assoc      1420      178.095
assoc      1430      182.049
assoc      1440      177.725
assoc      1450      179.253
assoc      1460      181.2
assoc      1470      182.122
assoc      1480      180.92
assoc      1490      179.435
assoc      1500      180.661
assoc      1510      180.382
assoc      1520      181.768
assoc      1530      180.342
assoc      1540      181.648
assoc      1550      181.178
assoc      1560      183.903
assoc      1570      182.538
assoc      1580      181.243
assoc      1590      183.426
assoc      1600      186.311
assoc      1610      184.986
assoc      1620      184.56
assoc      1630      183.473
assoc      1640      182.797
assoc      1650      183.016
assoc      1660      182.876
assoc      1670      183.612
assoc      1680      187.733
assoc      1690      184.976
assoc      1700      186.317
assoc      1710      186.696
assoc      1720      185.016
assoc      1730      191.803
assoc      1740      190.652
assoc      1750      186.834
assoc      1760      186.056
assoc      1770      185.311
assoc      1780      185.72
assoc      1790      187.316
assoc      1800      187.126
assoc      1810      188.341
assoc      1820      187.734
assoc      1830      187.41
assoc      1840      191.155
assoc      1850      190.326
assoc      1860      187.938
assoc      1870      188.124
assoc      1880      190.667
assoc      1890      188.345
assoc      1900      191.922
assoc      1910      188.809
assoc      1920      190.362
assoc      1930      192.464
assoc      1940      190.092
assoc      1950      189.901
assoc      1960      191.202
assoc      1970      192.834
assoc      1980      194.99
assoc      1990      189.845
assoc      2000      192.008
assoc      2010      190.895
assoc      2020      190.642
assoc      2030      190.83
assoc      2040      190.837
assoc      2050      191.898
assoc      2060      192.644
assoc      2070      193.474
assoc      2080      194.946
assoc      2090      193.121
assoc      2100      192.384
assoc      2110      193.209
assoc      2120      192.644
assoc      2130      194.561
assoc      2140      195.312
assoc      2150      196.94
assoc      2160      195.35
assoc      2170      192.86
assoc      2180      194.497
assoc      2190      195.129
assoc      2200      197.658
assoc      2210      195.149
assoc      2220      196.149
assoc      2230      194.4
assoc      2240      195.635
assoc      2250      197.024
assoc      2260      195.287
assoc      2270      199.246
assoc      2280      197.62
assoc      2290      195.04
assoc      2300      195.628
assoc      2310      196.539
assoc      2320      196.461
assoc      2330      196.582
assoc      2340      198.768
assoc      2350      195.822
assoc      2360      196.264
assoc      2370      199.489
assoc      2380      197.161
assoc      2390      198.202
assoc      2400      200.118
assoc      2410      198.228
assoc      2420      201.065
assoc      2430      199.882
assoc      2440      197.442
assoc      2450      197.082
assoc      2460      198.814
assoc      2470      197.586
assoc      2480      203.957
assoc      2490      199.327
assoc      2500      198.97
assoc      2510      199.327
assoc      2520      198.851
assoc      2530      198.456
assoc      2540      203.135
assoc      2550      207.656
assoc      2560      200.119
assoc      2570      203.257
assoc      2580      200.636
assoc      2590      200.199
assoc      2600      200.119
assoc      2610      199.563
assoc      2620      202.982
assoc      2630      200.596
assoc      2640      200.278
assoc      2650      199.722
assoc      2660      201.314
assoc      2670      198.886
assoc      2680      202.787
assoc      2690      201.193
assoc      2700      202.31
assoc      2710      203.701
assoc      2720      200.636
assoc      2730      205.763
assoc      2740      203.26
assoc      2750      203.221
assoc      2760      201.073
assoc      2770      204.138
assoc      2780      202.589
assoc      2790      200.517
assoc      2800      201.475
assoc      2810      199.283
assoc      2820      202.432
assoc      2830      206.978
assoc      2840      201.433
assoc      2850      202.466
assoc      2860      207.895
assoc      2870      204.823
assoc      2880      202.273
assoc      2890      203.708
assoc      2900      202.911
assoc      2910      205.425
assoc      2920      202.791
assoc      2930      204.143
assoc      2940      204.625
assoc      2950      203.911
assoc      2960      207.365
assoc      2970      204.982
assoc      2980      203.947
assoc      2990      202.154
assoc      3000      203.149
assoc      3010      204.629
assoc      3020      206.223
assoc      3030      203.992
assoc      3040      202.43
assoc      3050      204.79
assoc      3060      202.714
assoc      3070      205.465
assoc      3080      204.703
assoc      3090      205.502
assoc      3100      203.349
assoc      3110      204.622
assoc      3120      203.866
assoc      3130      204.109
assoc      3140      204.266
assoc      3150      202.595
assoc      3160      204.109
assoc      3170      203.831
assoc      3180      204.663
assoc      3190      205.022
assoc      3200      205.299
assoc      3210      205.54
assoc      3220      205.104
assoc      3230      204.948
assoc      3240      207.339
assoc      3250      206.302
assoc      3260      206.741
assoc      3270      206.975
assoc      3280      205.791
assoc      3290      207.228
assoc      3300      204.544
assoc      3310      207.057
assoc      3320      206.14
assoc      3330      204.663
assoc      3340      206.497
assoc      3350      204.358
assoc      3360      206.434
assoc      3370      205.762
assoc      3380      208.538
assoc      3390      205.505
assoc      3400      208.763
assoc      3410      206.352
assoc      3420      205.359
assoc      3430      206.984
assoc      3440      206.163
assoc      3450      207.282
assoc      3460      204.976
assoc      3470      206.476
assoc      3480      205.761
assoc      3490      207.441
assoc      3500      206.17
assoc      3510      208.038
assoc      3520      205.412
assoc      3530      208.194
assoc      3540      207.723
assoc      3550      207.402
assoc      3560      205.615
assoc      3570      206.656
assoc      3580      205.347
assoc      3590      206.807
assoc      3600      206.494
assoc      3610      213.426
assoc      3620      205.255
assoc      3630      207.137
assoc      3640      206.741
assoc      3650      208.45
assoc      3660      207.41
assoc      3670      207.3
assoc      3680      207.146
assoc      3690      206.587
assoc      3700      208.227
assoc      3710      206.385
assoc      3720      207.06
assoc      3730      206.225
assoc      3740      206.621
assoc      3750      206.576
assoc      3760      209.051
assoc      3770      208.892
assoc      3780      209.497
assoc      3790      208.011
assoc      3800      207.49
assoc      3810      209.773
assoc      3820      207.094
assoc      3830      207.766
assoc      3840      210.008
assoc      3850      210.606
assoc      3860      208.895
assoc      3870      209.015
assoc      3880      205.706
assoc      3890      207.018
assoc      3900      209.651
assoc      3910      209.268
assoc      3920      208.048
assoc      3930      208.37
assoc      3940      206.651
assoc      3950      207.903
assoc      3960      207.686
assoc      3970      208.367
assoc      3980      209.011
assoc      3990      209.326
assoc      4000      210.251
disassoc      0         294.573
disassoc      10        294.493
disassoc      20        294.02
disassoc      30        292.277
disassoc      40        293.548
disassoc      50        293.499
disassoc      60        294.596
disassoc      70        292.038
disassoc      80        290.898
disassoc      90        294.043
disassoc      100       290.694
disassoc      110       289.905
disassoc      120       290.232
disassoc      130       288.845
disassoc      140       292.18
disassoc      150       287.869
disassoc      160       288.14
disassoc      170       288.595
disassoc      180       286.963
disassoc      190       288.806
disassoc      200       288.136
disassoc      210       284.148
disassoc      220       283.556
disassoc      230       281.607
disassoc      240       283.846
disassoc      250       283.918
disassoc      260       283.918
disassoc      270       281.578
disassoc      280       279.354
disassoc      290       280
disassoc      300       282.446
disassoc      310       280.569
disassoc      320       281.725
disassoc      330       280.299
disassoc      340       277.83
disassoc      350       278.215
disassoc      360       279.969
disassoc      370       283.124
disassoc      380       279.283
disassoc      390       276.784
disassoc      400       278.312
disassoc      410       276.025
disassoc      420       277.129
disassoc      430       275.226
disassoc      440       275.868
disassoc      450       275.819
disassoc      460       272.86
disassoc      470       277.371
disassoc      480       272.452
disassoc      490       277.257
disassoc      500       275.187
disassoc      510       272.003
disassoc      520       272.566
disassoc      530       269.574
disassoc      540       273.028
disassoc      550       271.457
disassoc      560       272.942
disassoc      570       268.903
disassoc      580       269.98
disassoc      590       271.164
disassoc      600       268.664
disassoc      610       269.11
disassoc      620       269.173
disassoc      630       267.64
disassoc      640       267.021
disassoc      650       270.004
disassoc      660       268.309
disassoc      670       265.708
disassoc      680       270.484
disassoc      690       264.564
disassoc      700       266.759
disassoc      710       272.652
disassoc      720       266.43
disassoc      730       267.141
disassoc      740       264.43
disassoc      750       265.113
disassoc      760       264.283
disassoc      770       263.316
disassoc      780       263.083
disassoc      790       263.636
disassoc      800       264.531
disassoc      810       262.095
disassoc      820       261.797
disassoc      830       259.684
disassoc      840       259.945
disassoc      850       262.426
disassoc      860       259.311
disassoc      870       257.905
disassoc      880       258.017
disassoc      890       260.158
disassoc      900       260.008
disassoc      910       258.225
disassoc      920       258.294
disassoc      930       257.075
disassoc      940       258.879
disassoc      950       256.174
disassoc      960       256.748
disassoc      970       255.512
disassoc      980       255.88
disassoc      990       255.665
disassoc      1000      258.007
disassoc      1010      255.424
disassoc      1020      256.788
disassoc      1030      256.719
disassoc      1040      253.855
disassoc      1050      255.31
disassoc      1060      254.427
disassoc      1070      255.867
disassoc      1080      252.867
disassoc      1090      251.798
disassoc      1100      251.303
disassoc      1110      254.721
disassoc      1120      252.453
disassoc      1130      252.552
disassoc      1140      250.534
disassoc      1150      253.081
disassoc      1160      252.27
disassoc      1170      253.557
disassoc      1180      253.218
disassoc      1190      253.218
disassoc      1200      249.487
disassoc      1210      253.515
disassoc      1220      252.64
disassoc      1230      251.323
disassoc      1240      248.399
disassoc      1250      249.426
disassoc      1260      250.771
disassoc      1270      250.237
disassoc      1280      250.869
disassoc      1290      248.837
disassoc      1300      248.276
disassoc      1310      249.686
disassoc      1320      247.147
disassoc      1330      249.094
disassoc      1340      246.809
disassoc      1350      247.217
disassoc      1360      247.31
disassoc      1370      245.71
disassoc      1380      247.845
disassoc      1390      245.71
disassoc      1400      246.066
disassoc      1410      247.843
disassoc      1420      245.965
disassoc      1430      247.129
disassoc      1440      248.082
disassoc      1450      248.003
disassoc      1460      244.48
disassoc      1470      245.429
disassoc      1480      249.723
disassoc      1490      246.455
disassoc      1500      246.26
disassoc      1510      242.557
disassoc      1520      242.501
disassoc      1530      244.264
disassoc      1540      245.368
disassoc      1550      244.114
disassoc      1560      242.863
disassoc      1570      247.466
disassoc      1580      244.123
disassoc      1590      242.931
disassoc      1600      242.535
disassoc      1610      243.361
disassoc      1620      244.097
disassoc      1630      241.035
disassoc      1640      240.285
disassoc      1650      242.682
disassoc      1660      242.178
disassoc      1670      243.743
disassoc      1680      240.988
disassoc      1690      241.551
disassoc      1700      238.924
disassoc      1710      242.241
disassoc      1720      241.97
disassoc      1730      237.89
disassoc      1740      241.404
disassoc      1750      239.597
disassoc      1760      241.304
disassoc      1770      237.643
disassoc      1780      239.636
disassoc      1790      239.628
disassoc      1800      238.005
disassoc      1810      240.055
disassoc      1820      240.119
disassoc      1830      238.033
disassoc      1840      239.352
disassoc      1850      241.36
disassoc      1860      237.737
disassoc      1870      237.999
disassoc      1880      237.811
disassoc      1890      239.439
disassoc      1900      236.611
disassoc      1910      238.871
disassoc      1920      237.332
disassoc      1930      235.57
disassoc      1940      235.869
disassoc      1950      235.918
disassoc      1960      235.243
disassoc      1970      235.836
disassoc      1980      236.169
disassoc      1990      238.221
disassoc      2000      235.152
disassoc      2010      235.836
disassoc      2020      234.903
disassoc      2030      235.152
disassoc      2040      237.717
disassoc      2050      235.289
disassoc      2060      234.673
disassoc      2070      237.954
disassoc      2080      234.874
disassoc      2090      234.254
disassoc      2100      234.597
disassoc      2110      234.741
disassoc      2120      234.913
disassoc      2130      234.242
disassoc      2140      232.517
disassoc      2150      234.532
disassoc      2160      234.255
disassoc      2170      235.336
disassoc      2180      231.75
disassoc      2190      236.371
disassoc      2200      234.349
disassoc      2210      235.338
disassoc      2220      233.016
disassoc      2230      233.941
disassoc      2240      234.423
disassoc      2250      234.939
disassoc      2260      233.162
disassoc      2270      232.648
disassoc      2280      233.544
disassoc      2290      230.988
disassoc      2300      232.527
disassoc      2310      232.317
disassoc      2320      235.166
disassoc      2330      234.453
disassoc      2340      232.055
disassoc      2350      230.198
disassoc      2360      232.57
disassoc      2370      234.599
disassoc      2380      235.18
disassoc      2390      231.29
disassoc      2400      233.175
disassoc      2410      230.946
disassoc      2420      234.086
disassoc      2430      234.007
disassoc      2440      231.683
disassoc      2450      230.653
disassoc      2460      231.389
disassoc      2470      229.149
disassoc      2480      229.109
disassoc      2490      228.73
disassoc      2500      229.027
disassoc      2510      229.712
disassoc      2520      229.704
disassoc      2530      229.854
disassoc      2540      231.292
disassoc      2550      229.905
disassoc      2560      230.435
disassoc      2570      230.684
disassoc      2580      234.375
disassoc      2590      228.956
disassoc      2600      229.77
disassoc      2610      229.053
disassoc      2620      230.206
disassoc      2630      229.798
disassoc      2640      229.703
disassoc      2650      229.465
disassoc      2660      228.113
disassoc      2670      229.647
disassoc      2680      229.041
disassoc      2690      227.615
disassoc      2700      228.147
disassoc      2710      230.599
disassoc      2720      226.574
disassoc      2730      229.274
disassoc      2740      230.412
disassoc      2750      230.916
disassoc      2760      228.022
disassoc      2770      228.271
disassoc      2780      229.211
disassoc      2790      226.585
disassoc      2800      227.348
disassoc      2810      228.781
disassoc      2820      227.118
disassoc      2830      228.639
disassoc      2840      228.922
disassoc      2850      227.428
disassoc      2860      227.431
disassoc      2870      226.247
disassoc      2880      225.791
disassoc      2890      227.01
disassoc      2900      228.577
disassoc      2910      226.12
disassoc      2920      226.437
disassoc      2930      226.416
disassoc      2940      226.931
disassoc      2950      226.318
disassoc      2960      228
disassoc      2970      225.416
disassoc      2980      228.532
disassoc      2990      227.273
disassoc      3000      226.664
disassoc      3010      227.021
disassoc      3020      227.24
disassoc      3030      225.942
disassoc      3040      228.962
disassoc      3050      227.208
disassoc      3060      226.894
disassoc      3070      226.347
disassoc      3080      231.761
disassoc      3090      225.01
disassoc      3100      225.149
disassoc      3110      224.465
disassoc      3120      224.395
disassoc      3130      225.119
disassoc      3140      224.245
disassoc      3150      225.139
disassoc      3160      226.45
disassoc      3170      223.879
disassoc      3180      226.801
disassoc      3190      228
disassoc      3200      227.208
disassoc      3210      227.377
disassoc      3220      224.614
disassoc      3230      225.05
disassoc      3240      224.911
disassoc      3250      224.861
disassoc      3260      225.554
disassoc      3270      224.038
disassoc      3280      225.813
disassoc      3290      222.76
disassoc      3300      223.595
disassoc      3310      224.614
disassoc      3320      225.674
disassoc      3330      224.693
disassoc      3340      224.098
disassoc      3350      224.405
disassoc      3360      226.278
disassoc      3370      223.266
disassoc      3380      223.233
disassoc      3390      226.886
disassoc      3400      224.881
disassoc      3410      224.027
disassoc      3420      223.007
disassoc      3430      223.728
disassoc      3440      221.477
disassoc      3450      225.762
disassoc      3460      222.306
disassoc      3470      226.685
disassoc      3480      225.04
disassoc      3490      222.495
disassoc      3500      224.574
disassoc      3510      226.111
disassoc      3520      223.94
disassoc      3530      225.208
disassoc      3540      222.606
disassoc      3550      222.558
disassoc      3560      223.046
disassoc      3570      226.825
disassoc      3580      225.357
disassoc      3590      223.732
disassoc      3600      224.157
disassoc      3610      225.268
disassoc      3620      223.662
disassoc      3630      221.086
disassoc      3640      221.737
disassoc      3650      222.477
disassoc      3660      223.702
disassoc      3670      222.451
disassoc      3680      223.422
disassoc      3690      222.61
disassoc      3700      221.094
disassoc      3710      221.429
disassoc      3720      226.071
disassoc      3730      222.827
disassoc      3740      221.318
disassoc      3750      220.159
disassoc      3760      221.579
disassoc      3770      222.408
disassoc      3780      223.541
disassoc      3790      221.57
disassoc      3800      220.946
disassoc      3810      224.553
disassoc      3820      221.851
disassoc      3830      221.264
disassoc      3840      221.207
disassoc      3850      223.522
disassoc      3860      221.599
disassoc      3870      221.318
disassoc      3880      219.778
disassoc      3890      221.9
disassoc      3900      220.707
disassoc      3910      219.992
disassoc      3920      222.421
disassoc      3930      225.099
disassoc      3940      221.584
disassoc      3950      220.853
disassoc      3960      220.499
disassoc      3970      219.27
disassoc      3980      222.949
disassoc      3990      221.927
disassoc      4000      221.998
;

title1 'Systems of Differential Equations Example';

proc sort data=fit;
   by type time;
run;

%let k1f = 6.85e6  ;
%let k1r = 3.43e-4 ;
%let k2f = 1.8e7  ;
%let k2r = 2.1e-2 ;

%let qf  = 2.1e8  ;
%let qb  = 4.0e9  ;

%let dt  = 5.0e-7 ;
%let et  = 5.0e-8 ;
%let it  = 8.05e-6 ;

proc model data=fit;

   parameters qf  = 2.1e8
              qb  = 4.0e9
              k2f = 1.8e5
              k2r = 2.1e-3
              l   = 0;

              k1f = 6.85e6;
              k1r = 3.43e-4;

      /* Initial values for concentrations */
   control dt 5.0e-7
           et 5.0e-8
           it 8.05e-6;

      /* Association initial values --------------*/
   if type = 'assoc' and time=0 then do;
      ed = 0;
      /* solve quadratic equation ----------*/
      a = 1;
      b = -(&it+&et+(k2r/k2f));
      c = &it*&et;
      ei = (-b-(((b**2)-(4*a*c))**.5))/(2*a);
      d = &dt-ed;
      i = &it-ei;
      e = &et-ed-ei;
      end;

      /* Disassociation initial values ----------*/
   if type = 'disassoc' and time=0 then do;
         ei = 0;
         a = 1;
         b = -(&dt+&et+(&k1r/&k1f));
         c = &dt*&et;
         ed = (-b-(((b**2)-(4*a*c))**.5))/(2*a);
         d = &dt-ed;
         i = &it-ei;
         e = &et-ed-ei;
      end;

   if time ne 0 then do;

      dert.d = k1r* ed  - k1f *e *d;

      dert.ed = k1f* e *d - k1r*ed;

      dert.e = k1r* ed - k1f* e * d  + k2r * ei - k2f * e *i;

      dert.ei = k2f* e *i - k2r * ei;

      dert.i = k2r * ei - k2f* e *i;

      end;

   /* L - offset between curves  */
   if type = 'disassoc' then
      F = (qf*(d-ed)) + (qb*ed) -L;
   else
      F = (qf*(d-ed)) + (qb*ed);

   fit F / method=marquardt;
run;