From 8715bf4d29c342da859274ff0f270a52f3d51be9 Mon Sep 17 00:00:00 2001 From: Dorival Pedroso Date: Thu, 14 Mar 2024 20:44:34 +1000 Subject: [PATCH] [ode] Add one-transistor amplifier example --- russell_ode/README.md | 89 +- .../data/figures/amplifier1t_radau5.svg | 2991 +++++++++++++++++ .../reference/amplifier1t_mathematica.json | 1220 +++++++ russell_ode/examples/amplifier1t_radau5.rs | 111 + russell_ode/src/samples.rs | 87 +- 5 files changed, 4483 insertions(+), 15 deletions(-) create mode 100644 russell_ode/data/figures/amplifier1t_radau5.svg create mode 100644 russell_ode/data/reference/amplifier1t_mathematica.json create mode 100644 russell_ode/examples/amplifier1t_radau5.rs diff --git a/russell_ode/README.md b/russell_ode/README.md index a8fe522b..992ece61 100644 --- a/russell_ode/README.md +++ b/russell_ode/README.md @@ -660,4 +660,91 @@ The results are show below: ### One-transistor amplifier -TODO +This example corresponds to Fig 1.3 on page 377 and Fig 1.4 on page 379 of Reference #2. See also Eq (1.14) on page 377 of Reference #2. + +This is a differential-algebraic problem modelling the nodal voltages of a one-transistor amplifier. + +The DAE is expressed in the so-called *mass-matrix* form (ndim = 5): + +```text +M y' = f(x, y) + +with: y0(0)=0, y1(0)=Ub/2, y2(0)=Ub/2, y3(0)=Ub, y4(0)=0 +``` + +where the elements of the right-hand side function are: + +```text +f0 = (y0 - ue) / R +f1 = (2 y1 - UB) / S + γ g12 +f2 = y2 / S - g12 +f3 = (y3 - UB) / S + α g12 +f4 = y4 / S + +with: + +ue = A sin(ω x) +g12 = β (exp((y1 - y2) / UF) - 1) +``` + +Compared to Eq (1.14), we set all resistances Rᵢ to S, except the first one (R := R₀). + +The mass matrix is: + +```text + ┌ ┐ + │ -C1 C1 │ + │ C1 -C1 │ +M = │ -C2 │ + │ -C3 C3 │ + │ C3 -C3 │ + └ ┘ +``` + +and the Jacobian matrix is: + +```text + ┌ ┐ + │ 1/R │ + │ 2/S + γ h12 -γ h12 │ +J = │ -h12 1/S + h12 │ + │ α h12 -α h12 │ + │ 1/S │ + │ 1/S │ + └ ┘ + +with: + +h12 = β exp((y1 - y2) / UF) / UF +``` + +**Note:** In this library, only **Radau5** can solve such DAE. + +See the code [amplifier1t_radau5.rs](https://github.com/cpmech/russell/tree/main/russell_ode/examples/amplifier1t_radau5.rs) + +The output is given below: + +```text +y_russell = [-0.022267, 3.068709, 2.898349, 1.499405, -1.735090] +y_mathematica = [-0.022267, 3.068709, 2.898349, 1.499439, -1.735057] +Radau5: Radau method (Radau IIA) (implicit, order 5, embedded) +Number of function evaluations = 6007 +Number of Jacobian evaluations = 480 +Number of factorizations = 666 +Number of lin sys solutions = 1840 +Number of performed steps = 666 +Number of accepted steps = 481 +Number of rejected steps = 39 +Number of iterations (maximum) = 6 +Number of iterations (last step) = 1 +Last accepted/suggested stepsize = 0.00007705697843645314 +Max time spent on a step = 55.281µs +Max time spent on the Jacobian = 729ns +Max time spent on factorization = 249.11µs +Max time spent on lin solution = 241.201µs +Total time = 97.951021ms +``` + +The results are plotted below: + +![One-transistor Amplifier - Radau5](data/figures/amplifier1t_radau5.svg) diff --git a/russell_ode/data/figures/amplifier1t_radau5.svg b/russell_ode/data/figures/amplifier1t_radau5.svg new file mode 100644 index 00000000..55896ca3 --- /dev/null +++ b/russell_ode/data/figures/amplifier1t_radau5.svg @@ -0,0 +1,2991 @@ + + + + + + + + 2024-03-14T20:42:58.659377 + image/svg+xml + + + Matplotlib v3.6.3, https://matplotlib.org/ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/russell_ode/data/reference/amplifier1t_mathematica.json b/russell_ode/data/reference/amplifier1t_mathematica.json new file mode 100644 index 00000000..392519d5 --- /dev/null +++ b/russell_ode/data/reference/amplifier1t_mathematica.json @@ -0,0 +1,1220 @@ +{ + "x":[ + 0.0, + 1.0e-3, + 2.0e-3, + 3.0e-3, + 4.0e-3, + 5.0e-3, + 6.0e-3, + 7.0e-3, + 8.0e-3, + 9.000000000000001e-3, + 1.0e-2, + 1.1e-2, + 1.2e-2, + 1.3000000000000001e-2, + 1.4e-2, + 1.5e-2, + 1.6e-2, + 1.7e-2, + 1.8000000000000002e-2, + 1.9e-2, + 2.0e-2, + 2.1e-2, + 2.2e-2, + 2.3e-2, + 2.4e-2, + 2.5e-2, + 2.6000000000000002e-2, + 2.7e-2, + 2.8e-2, + 2.9e-2, + 3.0e-2, + 3.1e-2, + 3.2e-2, + 3.3e-2, + 3.4e-2, + 3.5e-2, + 3.6000000000000004e-2, + 3.7e-2, + 3.8e-2, + 3.9e-2, + 4.0e-2, + 4.1e-2, + 4.2e-2, + 4.3000000000000003e-2, + 4.4e-2, + 4.5e-2, + 4.6e-2, + 4.7e-2, + 4.8e-2, + 4.9e-2, + 5.0e-2, + 5.1000000000000004e-2, + 5.2000000000000005e-2, + 5.3e-2, + 5.4e-2, + 5.5e-2, + 5.6e-2, + 5.7e-2, + 5.8e-2, + 5.9000000000000004e-2, + 6.0e-2, + 6.1e-2, + 6.2e-2, + 6.3e-2, + 6.4e-2, + 6.5e-2, + 6.6e-2, + 6.7e-2, + 6.8e-2, + 6.9e-2, + 7.0e-2, + 7.100000000000001e-2, + 7.200000000000001e-2, + 7.3e-2, + 7.4e-2, + 7.5e-2, + 7.6e-2, + 7.7e-2, + 7.8e-2, + 7.9e-2, + 8.0e-2, + 8.1e-2, + 8.2e-2, + 8.3e-2, + 8.4e-2, + 8.5e-2, + 8.600000000000001e-2, + 8.700000000000001e-2, + 8.8e-2, + 8.9e-2, + 9.0e-2, + 9.1e-2, + 9.2e-2, + 9.3e-2, + 9.4e-2, + 9.5e-2, + 9.6e-2, + 9.7e-2, + 9.8e-2, + 9.9e-2, + 0.1, + 0.101, + 0.10200000000000001, + 0.10300000000000001, + 0.10400000000000001, + 0.105, + 0.106, + 0.107, + 0.108, + 0.109, + 0.11, + 0.111, + 0.112, + 0.113, + 0.114, + 0.115, + 0.116, + 0.117, + 0.11800000000000001, + 0.11900000000000001, + 0.12, + 0.121, + 0.122, + 0.123, + 0.124, + 0.125, + 0.126, + 0.127, + 0.128, + 0.129, + 0.13, + 0.131, + 0.132, + 0.133, + 0.134, + 0.135, + 0.136, + 0.137, + 0.138, + 0.139, + 0.14, + 0.14100000000000001, + 0.14200000000000002, + 0.14300000000000002, + 0.14400000000000002, + 0.145, + 0.146, + 0.147, + 0.148, + 0.149, + 0.15, + 0.151, + 0.152, + 0.153, + 0.154, + 0.155, + 0.156, + 0.157, + 0.158, + 0.159, + 0.16, + 0.161, + 0.162, + 0.163, + 0.164, + 0.165, + 0.166, + 0.167, + 0.168, + 0.169, + 0.17, + 0.171, + 0.17200000000000001, + 0.17300000000000001, + 0.17400000000000002, + 0.17500000000000002, + 0.176, + 0.177, + 0.178, + 0.179, + 0.18, + 0.181, + 0.182, + 0.183, + 0.184, + 0.185, + 0.186, + 0.187, + 0.188, + 0.189, + 0.19, + 0.191, + 0.192, + 0.193, + 0.194, + 0.195, + 0.196, + 0.197, + 0.198, + 0.199, + 0.2 + ], + "y0":[ + 0.0, + 0.19163123321667005, + 0.32185367911941054, + 0.3335960604269822, + 0.22156356989328296, + 2.8220796689663288e-2, + -0.17260350709122416, + -0.3052891386539197, + -0.3208689990238207, + -0.21128708301590177, + -1.9646365192382605e-2, + 0.1805568377307299, + 0.31261843258725863, + 0.3258768724679843, + 0.21510308403124437, + 2.281612798097559e-2, + -0.17710793265806787, + -0.3089262713111425, + -0.3239999935693118, + -0.21390588121750806, + -2.183390978444637e-2, + 0.17872962198567616, + 0.3110919993525171, + 0.3246008834677472, + 0.2140349944794215, + 2.1922590884010797e-2, + -0.17785273664260093, + -0.3095276892691252, + -0.3245174973765986, + -0.21433880893641322, + -2.219554644903774e-2, + 0.17842755522289921, + 0.3108396554171122, + 0.3243900605329387, + 0.21385863130346497, + 2.1775038086288823e-2, + -0.17797573192859859, + -0.3096270054719691, + -0.32460295275303813, + -0.21441030318818705, + -2.225526702199502e-2, + 0.17837766853446904, + 0.31079798125737845, + 0.3243552252705964, + 0.21382947415638637, + 2.1750644693117698e-2, + -0.1779960653623523, + -0.3096434251784288, + -0.32461707984533517, + -0.21442212082478654, + -2.2265138739697748e-2, + 0.17836942449407503, + 0.3107910941805478, + 0.32434946933218667, + 0.21382465623799582, + 2.174661397691869e-2, + -0.1779994253672666, + -0.3096461384007443, + -0.32461941411601647, + -0.21442407404728434, + -2.2266770123088607e-2, + 0.17836806269927047, + 0.31078995374478746, + 0.32434851699480577, + 0.21382385941913695, + 2.174594910244149e-2, + -0.17799997968327477, + -0.30964658577676113, + -0.32461963510442166, + -0.2144239795092126, + -2.226641282093319e-2, + 0.1783683495652274, + 0.3107901983730448, + 0.3243487203513457, + 0.21382402937304035, + 2.1746089601628542e-2, + -0.17799986228697626, + -0.30964649138636213, + -0.3246197182253148, + -0.21442432875230935, + -2.2266982460730687e-2, + 0.17836788520443284, + 0.3107898079103974, + 0.32434825990459754, + 0.21382354655611355, + 2.1745697685687643e-2, + -0.1780001884785888, + -0.30964676101942845, + -0.3246199425596532, + -0.21442451607731042, + -2.2267139406040557e-2, + 0.17836775376684744, + 0.31078969830768266, + 0.3243483023585903, + 0.21382367971301652, + 2.1745797027226375e-2, + -0.17800010603672206, + -0.3096466881064059, + -0.32461988713464335, + -0.21442446998461506, + -2.226710115357194e-2, + 0.17836778490305594, + 0.31078965739623077, + 0.3243481383032177, + 0.2138234615765023, + 2.1745627335258876e-2, + -0.17800024711591345, + -0.30964680878481887, + -0.3246199832951104, + -0.21442455099173957, + -2.226716863142414e-2, + 0.1783677301951348, + 0.3107896760827645, + 0.3243480770428223, + 0.2138234100428865, + 2.1745584187124047e-2, + -0.17800028284498673, + -0.3096468375739029, + -0.3246200082462933, + -0.2144245708614127, + -2.226718531961889e-2, + 0.1783677163007287, + 0.3107896665725322, + 0.3243482754831915, + 0.21382365727366714, + 2.1745778262873647e-2, + -0.17800012200417772, + -0.3096467008529888, + -0.3246198981349176, + -0.21442447869217532, + -2.2267108221413416e-2, + 0.17836778065040942, + 0.3107897204018311, + 0.32434832086881754, + 0.21382369509113344, + 2.174581015266233e-2, + -0.17800009544588138, + -0.30964667950315944, + -0.32461987981426454, + -0.21442446415021096, + -2.2267097207446927e-2, + 0.17836778967245462, + 0.3107897280765495, + 0.3243481181774326, + 0.2138234443497565, + 2.174561313457391e-2, + -0.1780002588157916, + -0.30964681842426467, + -0.32461999166815314, + -0.21442455683650197, + -2.226717352074859e-2, + 0.17836772610714258, + 0.3107896746776564, + 0.3243482824265437, + 0.21382366307765027, + 2.174578871150656e-2, + -0.17800011200529972, + -0.30964669296675723, + -0.32461989115780415, + -0.21442447169081946, + -2.226710325271905e-2, + 0.17836778406399334, + 0.31078972216518685, + 0.3243483218240611, + 0.21382369516099045, + 2.174581023819502e-2, + -0.1780000953706733, + -0.3096466795362259, + -0.324619879701329, + -0.21442446346575783, + -2.2267095411230155e-2, + 0.17836779093851032, + 0.3107897290553566, + 0.3243483281848601, + 0.21382370216933316, + 2.174581736682977e-2, + -0.17800008974432593, + -0.3096466750377847, + -0.32461987574706114, + -0.21442446013579616, + -2.2267092718279556e-2, + 0.17836779064296604, + 0.31078972878311406, + 0.3243483259170929, + 0.21382369993780845, + 2.1745816419619774e-2, + -0.17800009012072013, + -0.3096466755584529, + -0.32461987593547964, + -0.21442446111595126, + -2.2267093340786855e-2, + 0.1783677925874619, + 0.31078973061269577, + 0.3243483291338532, + 0.21382370208184195, + 2.1745816943518494e-2, + -0.1780000899502744, + -0.30964667500175497, + -0.3246198759212826, + -0.21442445992979184, + -2.22670928295124e-2 + ], + "y1":[ + 3.0, + 3.168238648380042, + 3.2452403388220494, + 3.202121593005678, + 3.0585653225078024, + 2.872508886258942, + 2.7183210816701195, + 2.6567803578337776, + 2.7105244507232444, + 2.8630846620833936, + 3.0567809231266154, + 3.218527579121905, + 3.287248791312676, + 3.2371921918885067, + 3.0878242376215046, + 2.8969055203903507, + 2.7386747720960543, + 2.673920390313233, + 2.724802567671607, + 2.875006535128904, + 3.066737309978718, + 3.226842309778451, + 3.2941920179773247, + 3.2429888283328805, + 3.092660745528837, + 2.9009381558006773, + 2.742038874941794, + 2.6767527634976993, + 2.727163063095718, + 2.876977415947597, + 3.068383261600366, + 3.2282168609176636, + 3.295339840203077, + 3.243946899990163, + 3.0934594206648263, + 2.9016040681370674, + 2.742594385746785, + 2.67722045848238, + 2.72755286988428, + 2.877302883352526, + 3.0686550750306303, + 3.22844386981916, + 3.2955294059134768, + 3.2441051587225975, + 3.0935914476884374, + 2.9017141548001804, + 2.7426862212329692, + 2.677297774721051, + 2.7276173113576467, + 2.8773566851353007, + 3.0687000041933272, + 3.2284813815609237, + 3.2955607295488836, + 3.244131308885903, + 3.093613264238205, + 2.9017323452553803, + 2.742701395668836, + 2.6773105501726597, + 2.7276279596427893, + 2.877365575984673, + 3.0687074296353067, + 3.2284875801104724, + 3.2955659034997056, + 3.244135631297273, + 3.0936168727072326, + 2.901735345832867, + 2.742703899075404, + 2.677312657856236, + 2.7276295568828846, + 2.8773655944611654, + 3.068705816557579, + 3.2284862584347067, + 3.295564804771108, + 3.244134710904433, + 3.0936161025653317, + 2.901734711804185, + 2.742703370384386, + 2.677312212439086, + 2.7276293447928177, + 2.8773667340960394, + 3.068708397137011, + 3.228488388150614, + 3.2955665790956417, + 3.2441364344663817, + 3.0936182036594184, + 2.9017364702439266, + 2.742704836642212, + 2.6773134377862204, + 2.7276303697451065, + 2.8773675884807526, + 3.0687091104143107, + 3.228488985950396, + 3.295567080357178, + 3.244136608407684, + 3.0936176858941806, + 2.9017360325616637, + 2.7427044718975604, + 2.677313139999064, + 2.7276301184142726, + 2.8773673783460754, + 3.0687089360076536, + 3.2284888471994115, + 3.2955667707831404, + 3.244136863116259, + 3.093618583516717, + 2.901736787260218, + 2.74270510053401, + 2.6773136594349567, + 2.7276305545216206, + 2.877367745858061, + 3.0687092450957585, + 3.228489096563672, + 3.2955671729010527, + 3.2441371380328796, + 3.093618816537504, + 2.90173698183332, + 2.7427052633583484, + 2.6773137964217297, + 2.7276306685099776, + 2.8773678379319327, + 3.0687093184625334, + 3.228489156935781, + 3.29556722254832, + 3.244136730851819, + 3.093617788162887, + 2.9017361168693325, + 2.7427045419749043, + 2.6773131990988963, + 2.727630167410159, + 2.877367419099297, + 3.0687089687043922, + 3.2284888650015677, + 3.295566977360107, + 3.244136525468458, + 3.0936176153986974, + 2.9017359726901084, + 2.74270442275112, + 2.677313098706774, + 2.727630083695085, + 2.877367350661222, + 3.0687089179731304, + 3.228488822952391, + 3.295566945201563, + 3.244136944349928, + 3.093618660745421, + 2.901736851753002, + 2.742705155175675, + 2.677313705203849, + 2.727630592622157, + 2.8773677745955992, + 3.0687092656181276, + 3.228489113191071, + 3.2955671872490813, + 3.2441367003278883, + 3.0936177608359294, + 2.9017360630750337, + 2.7427044965104668, + 2.6773131602607134, + 2.7276301350826357, + 2.877367392769922, + 3.0687089458143593, + 3.228488845690917, + 3.29556696179161, + 3.2441365132629114, + 3.0936176147712477, + 2.9017359724666374, + 2.742704421736051, + 2.6773130977200363, + 2.7276300829226896, + 2.8773673486454565, + 3.0687089101343568, + 3.2284888169960118, + 3.295566938832484, + 3.24413649269672, + 3.0936175842357265, + 2.9017359418049535, + 2.7427043961569977, + 2.6773130762290553, + 2.7276300650624536, + 2.8773673336119745, + 3.0687088971322396, + 3.228488819734261, + 3.295566942444655, + 3.244136495416661, + 3.093617592537063, + 2.90173594390016, + 2.742704397647424, + 2.677313077116184, + 2.72763006591802, + 2.877367337291712, + 3.0687089008903286, + 3.2284888086926586, + 3.295566932415329, + 3.244136487012465, + 3.0936175845354406, + 2.9017359438135055, + 2.7427043972072984, + 2.6773130772430513, + 2.727630065805837, + 2.877367334589793, + 3.068708898325494 + ], + "y2":[ + 3.0, + 3.0014610192558364, + 3.088940548522638, + 3.065418109803608, + 2.9545637314628492, + 2.807779837225976, + 2.6598638305453495, + 2.5334672619469853, + 2.5500912345892788, + 2.6942525195893685, + 2.886310660291691, + 3.052183244655611, + 3.131593794589612, + 3.1015141756739895, + 2.985877675174497, + 2.8360795454185457, + 2.685907135057466, + 2.5548228971187257, + 2.5645981695144773, + 2.706294692543181, + 2.8963595094850496, + 3.0605868858409946, + 3.138645043179063, + 3.1074825034015476, + 2.99106124821594, + 2.840773064473547, + 2.69023678119389, + 2.5584121907786352, + 2.566997348771377, + 2.7082855214039294, + 2.8980207796669677, + 3.061976165580377, + 3.139810762656469, + 3.1084701051363046, + 2.9919184936871033, + 2.8415491550554037, + 2.69095291563894, + 2.559006820959752, + 2.5673935717967886, + 2.7086142824440547, + 2.8982951232433942, + 3.0622056070470043, + 3.1400032885437734, + 3.108633093841691, + 2.9920600639698027, + 2.8416773920585467, + 2.6910712634114065, + 2.559105142250299, + 2.567459074506787, + 2.708668630659485, + 2.8983404705395563, + 3.062243520684236, + 3.140035096977851, + 3.1086600278935097, + 2.992083458775125, + 2.841698584870206, + 2.691090819136992, + 2.5591213881953663, + 2.5674698984464155, + 2.7086776116946574, + 2.8983479651287127, + 3.0622497845829058, + 3.140040350978672, + 3.1086644771872396, + 2.992087333630618, + 2.8417020694977384, + 2.6910940383903106, + 2.5591240752095903, + 2.5674722319822285, + 2.7086780201740734, + 2.898346347636537, + 3.0622484331652937, + 3.1400392363015617, + 3.108663529514867, + 2.9920865010360584, + 2.8417013400115327, + 2.691093363917457, + 2.5591235026957024, + 2.5674713041867085, + 2.7086787800553993, + 2.898348942652621, + 3.062250601320543, + 3.1400410388563946, + 3.108664153587642, + 2.992087775157333, + 2.841702803566606, + 2.691094769810588, + 2.5591248010665706, + 2.567472347833715, + 2.7086796443793464, + 2.898349661875436, + 3.0622512060587184, + 3.140041549007966, + 3.1086654827902205, + 2.992088199277277, + 2.84170289430752, + 2.691094795859683, + 2.55912468916443, + 2.5674720943129747, + 2.7086794319249172, + 2.898349485579834, + 3.062251067617674, + 3.140040519377192, + 3.108664209854291, + 2.992088123209928, + 2.841703136938321, + 2.6910950795438775, + 2.5591250644913996, + 2.56747253457024, + 2.7086798021822616, + 2.8983497984030597, + 3.0622513172212162, + 3.140041630108178, + 3.10866448087176, + 2.992088367889477, + 2.8417033608394417, + 2.691095289049845, + 2.559125245571858, + 2.5674726506682033, + 2.7086798961457834, + 2.8983498722226857, + 3.0622513776328963, + 3.140041691213356, + 3.1086656076890997, + 2.9920883080173035, + 2.841702974448136, + 2.6910948729120263, + 2.559124759719089, + 2.567472142507831, + 2.708679473299202, + 2.898349518591936, + 3.062251084337811, + 3.1400414429888213, + 3.108665397256375, + 2.9920881204284777, + 2.8417028016033923, + 2.6910947280899604, + 2.5591246321694774, + 2.5674720572896423, + 2.708679402746572, + 2.8983494670524204, + 3.062251040879985, + 3.1400414135211108, + 3.108664261927427, + 2.992088193660677, + 2.84170320333418, + 2.6910951501914875, + 2.5591251247897784, + 2.5674725733132244, + 2.7086798322993744, + 2.898349818150705, + 3.0622513343666258, + 3.1400416563267073, + 3.108665578025622, + 2.992088277928836, + 2.8417028544801033, + 2.6910947658517768, + 2.559124700023755, + 2.567472109282841, + 2.7086794464220882, + 2.8983494954111397, + 3.0622510636252214, + 3.1400414261062717, + 3.108665386642899, + 2.9920881196486464, + 2.8417028044050294, + 2.691094714489794, + 2.559124624882496, + 2.5674720563797506, + 2.708679401939023, + 2.8983494593429153, + 3.062251034728043, + 3.140041403347032, + 3.1086653633528196, + 2.992088093503792, + 2.841702766302079, + 2.69109468135245, + 2.559124595691871, + 2.567472038398548, + 2.7086793867418657, + 2.8983494479247263, + 3.062251038854665, + 3.1400414076454344, + 3.1086653657586854, + 2.9920880989980256, + 2.8417027536351034, + 2.6910946668437283, + 2.5591245873976733, + 2.5674720393214927, + 2.7086793899225747, + 2.8983494501911378, + 3.062251026355273, + 3.140041397111356, + 3.108665357127274, + 2.9920880903690663, + 2.8417027716959757, + 2.6910946872095054, + 2.5591246034735264, + 2.5674720390378925, + 2.7086793879955153, + 2.898349447438977 + ], + "y3":[ + 6.0, + 3.25673471201135, + 4.117548762295695, + 5.057129946809468, + 5.661960512889709, + 5.851145305299047, + 5.863842393284494, + 5.393714712000808, + 3.749668695861498, + 2.891854773386496, + 2.64674589947657, + 3.0568172343525326, + 3.9225934420564617, + 4.855840326829122, + 5.4508435842170195, + 5.633464708449331, + 5.651222263719712, + 5.253922287382534, + 3.5569619920646987, + 2.6982245737482073, + 2.4546839258341455, + 2.866671222865716, + 3.7344747607016022, + 4.669486702134018, + 5.265652872030741, + 5.449941740562292, + 5.47118227424525, + 5.088530115790203, + 3.385669396278748, + 2.5292858567391927, + 2.2885162414315956, + 2.703286502917022, + 3.5738458920696448, + 4.511558893839967, + 5.110224188985182, + 5.297071404143296, + 5.321133982308646, + 4.943107682673369, + 3.2414450739906546, + 2.387575482120113, + 2.149351871214937, + 2.566631647507131, + 3.4396587076169456, + 4.379788247970368, + 4.980808865293592, + 5.169983160458308, + 5.196380427631806, + 4.82095221572049, + 3.121288558147732, + 2.2696019238547542, + 2.0335337395619173, + 2.4529312694518093, + 3.3280377025280976, + 4.2702080837448255, + 4.873228925589861, + 5.064369259564147, + 5.092704793601098, + 4.719229892996888, + 3.0213926686901287, + 2.1715352830755, + 1.9372649662756958, + 2.358427535615296, + 3.2352668860793212, + 4.179137971880075, + 4.783828384675777, + 4.976607700851542, + 5.0065535756054365, + 4.634667927599898, + 2.9384328114005043, + 2.0900831675988103, + 1.857263965849783, + 2.279891846549828, + 3.158174047675323, + 4.103460048075946, + 4.709539100871492, + 4.903681836660804, + 4.934965751795411, + 4.56439170881905, + 2.8693900541382207, + 2.022321204653687, + 1.7907887542869854, + 2.214638613201096, + 3.094116186362104, + 4.040541320332112, + 4.64780156321058, + 4.843083413015688, + 4.875479620396921, + 4.505998299631554, + 2.8120684406414536, + 1.9660509842842289, + 1.7355506193478092, + 2.1604140120016866, + 3.0408864127447237, + 3.9883237276262804, + 4.596514606623175, + 4.792730573407776, + 4.826050276489688, + 4.457478426517079, + 2.7644357083788007, + 1.9192921354112555, + 1.6896497254144296, + 2.1153555677205156, + 2.996605947183292, + 3.944855827418945, + 4.553882871190717, + 4.750888259704054, + 4.7849760769622245, + 4.417155952751167, + 2.724855693076282, + 1.880438317624667, + 1.651508960148814, + 2.077914222898954, + 2.9598994941874555, + 3.9087756622655587, + 4.5184649291077665, + 4.7161201151051415, + 4.750845907915054, + 4.383652661501687, + 2.691966581227511, + 1.8481526993431574, + 1.619815570870033, + 2.0468022862637527, + 2.929359199649866, + 3.8788429089631125, + 4.489042527144822, + 4.687230339126797, + 4.722485761243332, + 4.355814313560746, + 2.6646367198693195, + 1.8213243911317303, + 1.5934794789459117, + 2.020949715777339, + 2.903980897659619, + 3.853930289366357, + 4.464587027652546, + 4.663223638520934, + 4.698919571782832, + 4.332680101984467, + 2.6419272661408537, + 1.799031492389788, + 1.5715957458970533, + 1.9994674586561165, + 2.882893120730747, + 3.833180674038447, + 4.444257581060411, + 4.643274460152693, + 4.679336794799186, + 4.313455162657671, + 2.623057425332706, + 1.7805079780554494, + 1.5534119053521607, + 1.9816171890960395, + 2.865370097137491, + 3.8160278852384812, + 4.427379984617659, + 4.62669912262218, + 4.663065223265395, + 4.297483758011259, + 2.607376655425777, + 1.765114978329748, + 1.538301406935675, + 1.966783955178119, + 2.8508090229100618, + 3.80173396368164, + 4.413348536885071, + 4.612925287898606, + 4.649544073373872, + 4.284210452974201, + 2.594347037475048, + 1.7523244262039614, + 1.5257455565654554, + 1.954458491513565, + 2.8387097388290616, + 3.7898568238536647, + 4.4016891944647565, + 4.601479853935656, + 4.638308618549464, + 4.273181044846794, + 2.5835201104693106, + 1.7416961409585536, + 1.51531248116251, + 1.9442168021884725, + 2.828655870790848, + 3.7799872909192516, + 4.392000799433407, + 4.591969223357147, + 4.628972508959287, + 4.264016116335328, + 2.5745234703957074, + 1.7328645211308336, + 1.5066427823183761, + 1.9357061867319163, + 2.820301540048434, + 3.771786306249244, + 4.383950284838001, + 4.584066432174485, + 4.621214703280103, + 4.2564008873825285, + 2.567047712392508, + 1.72552602342074, + 1.4994388165471992 + ], + "y4":[ + 0.0, + -2.688717453566241, + -1.7448757696025976, + -0.7595061737742214, + -0.13954979725839428, + 5.0343984180890596e-2, + 6.067077242219998e-2, + -0.4072839950684525, + -2.0035631815070243, + -2.7712187089552143, + -2.90903469959208, + -2.398877913910467, + -1.4608517199613893, + -0.4920353588721492, + 0.1085394338107019, + 0.28299835731304424, + 0.2898784876165033, + -0.11445034512233938, + -1.7735054563915458, + -2.550406870782114, + -2.6947104121740244, + -2.1904647766558067, + -1.2580261884614234, + -0.2948540050354282, + 0.2996908582376115, + 0.4688415889970111, + 0.4723904796799939, + 7.588831964316665e-2, + -1.5959306795481907, + -2.376973034370672, + -2.5248620881467487, + -2.0240600293120417, + -1.094966544231726, + -0.1350710869109124, + 0.4561176429172534, + 0.6220930215475249, + 0.6228411473733575, + 0.2254236249245619, + -1.4506653409543047, + -2.2345208279272457, + -2.385090049258476, + -1.8869065466529016, + -0.9603775008385681, + -3.0034764259296073e-3, + 0.5856953973958074, + 0.7492439080348212, + 0.7476609222096344, + 0.34825774094956424, + -1.3303366586386858, + -2.1164241971914834, + -2.2691714717345315, + -1.7731238201661872, + -0.8486906494267208, + 0.1066258783832785, + 0.693302195871194, + 0.8548681604864885, + 0.8513474565490029, + 0.4500921319785882, + -1.230412248993869, + -2.018337285816265, + -2.1728860754108976, + -1.6786066659068923, + -0.7559089786042534, + 0.19770354205066898, + 0.782707251488033, + 0.9426313938544162, + 0.937500457551053, + 0.5346728103099889, + -1.1473329715505807, + -1.9367972963191906, + -2.0928861398687, + -1.6000770749852893, + -0.6788183320218945, + 0.27338016461535536, + 0.8569954667084551, + 1.01555693258891, + 1.0090879451604562, + 0.6049449291734457, + -1.0784063253817395, + -1.8691208608585312, + -2.026407460971851, + -1.534815963933954, + -0.614756623294818, + 0.3362298762686326, + 0.9187189389424409, + 1.0761543308035657, + 1.0685739292783667, + 0.6633407121940977, + -1.0210816690634246, + -1.8128483943828004, + -1.9711678871342118, + -1.4805899384529715, + -0.5615256545676945, + 0.3885199737730419, + 0.9700219212741125, + 1.1265089941653454, + 1.118004276395733, + 0.7118665966147016, + -0.9734492996655012, + -1.7660900664267143, + -1.9252674589356402, + -1.4355314317724805, + -0.5173429363684963, + 0.4318919916133658, + 1.0126370905954827, + 1.168349572538896, + 1.15907756475777, + 0.7521843957640715, + -0.9338686090728946, + -1.7272356429941473, + -1.8871257946127216, + -1.3980899999372227, + -0.480540285848618, + 0.4679719131688229, + 1.0480552309703923, + 1.2031178000322753, + 1.1932079064634837, + 0.7856891586791774, + -0.9009791561717678, + -1.6949496135650572, + -1.8554323083622777, + -1.3669780931571056, + -0.449998425176401, + 0.4980009679252728, + 1.0774941172494887, + 1.2320092076605142, + 1.2215687511006148, + 0.8135310485497868, + -0.8736504232062651, + -1.6681221895642053, + -1.8290971334510946, + -1.3411258008543654, + -0.42462054067745414, + 0.5229132632786169, + 1.1019493438825734, + 1.2560158155168688, + 1.2451349231788373, + 0.8366644121219937, + -0.850941212117499, + -1.6458297666310866, + -1.8072135962927074, + -1.319643832446047, + -0.403532400421603, + 0.5435653142144232, + 1.1222622640939655, + 1.2759634061777783, + 1.2647168993431581, + 0.8558856581610602, + -0.8320702650087038, + -1.6273049845573953, + -1.7890289245702986, + -1.3017928959771223, + -0.3860092791228095, + 0.5608160357432552, + 1.1391565916716557, + 1.2925401980997533, + 1.2809891602459054, + 0.8718609478194118, + -0.8163905231736043, + -1.6119128605824953, + -1.7739191337276439, + -1.2869603140594088, + -0.3714488402349656, + 0.5751093114969689, + 1.1531878275032492, + 1.3063141824205016, + 1.2945103778372358, + 0.8851338252667723, + -0.8033610246276064, + -1.5991224297653912, + -1.7613633616731068, + -1.274634904260414, + -0.3593495461161772, + 0.5869867513655714, + 1.1648472127080056, + 1.3177596172609607, + 1.3057457886853336, + 0.8961629857232228, + -0.7925341127121318, + -1.588494182373047, + -1.7509299205106372, + -1.2643929339641702, + -0.349295590794659, + 0.5968561007225635, + 1.1745355295162994, + 1.3272701667932219, + 1.3150818557692003, + 0.9053276366346846, + -0.7835374834152088, + -1.5796626746765365, + -1.7422605769968647, + -1.2558826495648876, + -0.3409413224220096, + 0.6050571445640937, + 1.1825860972287858, + 1.3351730420637973, + 1.322839713631253, + 0.9129433871775408, + -0.7760617267748756, + -1.5723239900112043, + -1.735056659198201 + ] +} \ No newline at end of file diff --git a/russell_ode/examples/amplifier1t_radau5.rs b/russell_ode/examples/amplifier1t_radau5.rs new file mode 100644 index 00000000..64bd3483 --- /dev/null +++ b/russell_ode/examples/amplifier1t_radau5.rs @@ -0,0 +1,111 @@ +use plotpy::{Curve, Plot}; +use russell_lab::StrError; +use russell_ode::prelude::*; +use serde::Deserialize; +use std::{env, fs::File, io::BufReader, path::Path}; + +// This example solves the DAE representing a one-transistor amplifier +// +// This example corresponds to Fig 1.3 on page 377 and Fig 1.4 on page 379 of the reference. +// See also Eq (1.14) on page 377 of the reference. +// +// # Reference +// +// * Hairer E, Wanner G (2002) Solving Ordinary Differential Equations II. +// Stiff and Differential-Algebraic Problems. Second Revised Edition. +// Corrected 2nd printing 2002. Springer Series in Computational Mathematics, 614p + +fn main() -> Result<(), StrError> { + // get the ODE system + let (system, mut data, mut args) = Samples::amplifier1t(); + + // solver + let params = Params::new(Method::Radau5); + let mut solver = OdeSolver::new(params, &system)?; + + // enable dense output + let mut out = Output::new(); + let h_out = 0.0001; + let selected_y_components = &[0, 4]; + out.enable_dense(h_out, selected_y_components)?; + + // solve the problem + let x1 = 0.2; + solver.solve(&mut data.y0, data.x0, x1, None, Some(&mut out), &mut args)?; + + // print the results and stats + let y_ref = &[ + -0.0222670928295124, + 3.068708898325494, + 2.898349447438977, + 1.4994388165471992, + -1.735056659198201, + ]; + println!("y_russell = {:.6?}", data.y0.as_data()); + println!("y_mathematica = {:.6?}", y_ref); + println!("{}", solver.bench()); + + // load reference data (from Mathematica) + let math = ReferenceData::read("data/reference/amplifier1t_mathematica.json")?; + + // plot the results + let mut curve1 = Curve::new(); + let mut curve2 = Curve::new(); + let mut curve3 = Curve::new(); + let mut curve4 = Curve::new(); + curve1.set_label("y0: russell"); + curve2.set_label("y0: mathematica"); + curve3.set_label("y4: russell"); + curve4.set_label("y4: mathematica"); + + let blue = "#307BC2"; + let red = "#C23048"; + + curve1.set_line_color(&blue); + curve2 + .set_marker_color(&blue) + .set_marker_line_color(&blue) + .set_marker_style(".") + .set_line_style("None"); + curve3.set_line_color(&red); + curve4 + .set_marker_color(&red) + .set_marker_line_color(&red) + .set_marker_style("+") + .set_line_style("None"); + + curve1.draw(&out.dense_x, out.dense_y.get(&0).unwrap()); + curve2.draw(&math.x, &math.y0); + curve3.draw(&out.dense_x, out.dense_y.get(&4).unwrap()); + curve4.draw(&math.x, &math.y4); + + // save figure + let mut plot = Plot::new(); + plot.add(&curve1) + .add(&curve2) + .add(&curve3) + .add(&curve4) + .legend() + .set_title("One-transistor Amplifier - Radau5") + .set_figure_size_points(800.0, 400.0) + .grid_and_labels("$y_0$", "$y_1$") + .save("/tmp/russell_ode/amplifier1t_radau5.svg") +} + +#[derive(Deserialize)] +struct ReferenceData { + pub x: Vec, + pub y0: Vec, + pub y4: Vec, +} + +impl ReferenceData { + pub fn read(rel_path: &str) -> Result { + let full_path = format!("{}/{}", env::var("CARGO_MANIFEST_DIR").unwrap(), rel_path); + let path = Path::new(full_path.as_str()).to_path_buf(); + let input = File::open(path).map_err(|_| "cannot open file")?; + let buffered = BufReader::new(input); + let data = serde_json::from_reader(buffered).map_err(|_| "cannot parse JSON file")?; + Ok(data) + } +} diff --git a/russell_ode/src/samples.rs b/russell_ode/src/samples.rs index dc0bebde..dc41caa7 100644 --- a/russell_ode/src/samples.rs +++ b/russell_ode/src/samples.rs @@ -714,6 +714,65 @@ impl Samples { /// Returns the one-transistor amplifier problem described by Hairer-Wanner, Part II, page 376 /// + /// This example corresponds to Fig 1.3 on page 377 and Fig 1.4 on page 379 of the reference. + /// See also Eq (1.14) on page 377 of the reference. + /// + /// This is a differential-algebraic problem modelling the nodal voltages of a one-transistor amplifier. + /// + /// The DAE is expressed in the so-called *mass-matrix* form (ndim = 5): + /// + /// ```text + /// M y' = f(x, y) + /// + /// with: y0(0)=0, y1(0)=Ub/2, y2(0)=Ub/2, y3(0)=Ub, y4(0)=0 + /// ``` + /// + /// where the elements of the right-hand side function are: + /// + /// ```text + /// f0 = (y0 - ue) / R + /// f1 = (2 y1 - UB) / S + γ g12 + /// f2 = y2 / S - g12 + /// f3 = (y3 - UB) / S + α g12 + /// f4 = y4 / S + /// + /// with: + /// + /// ue = A sin(ω x) + /// g12 = β (exp((y1 - y2) / UF) - 1) + /// ``` + /// + /// Compared to Eq (1.14), we set all resistances Rᵢ to S, except the first one (R := R₀). + /// + /// The mass matrix is: + /// + /// ```text + /// ┌ ┐ + /// │ -C1 C1 │ + /// │ C1 -C1 │ + /// M = │ -C2 │ + /// │ -C3 C3 │ + /// │ C3 -C3 │ + /// └ ┘ + /// ``` + /// + /// and the Jacobian matrix is: + /// + /// ```text + /// ┌ ┐ + /// │ 1/R │ + /// │ 2/S + γ h12 -γ h12 │ + /// J = │ -h12 1/S + h12 │ + /// │ α h12 -α h12 │ + /// │ 1/S │ + /// │ 1/S │ + /// └ ┘ + /// + /// with: + /// + /// h12 = β exp((y1 - y2) / UF) / UF + /// ``` + /// /// # Output /// /// Returns `(system, data, args)` where: @@ -742,9 +801,9 @@ impl Samples { // constants const ALPHA: f64 = 0.99; const GAMMA: f64 = 1.0 - ALPHA; - const C: f64 = 0.4; - const D: f64 = 200.0 * PI; const BETA: f64 = 1e-6; + const A: f64 = 0.4; + const OM: f64 = 200.0 * PI; const UB: f64 = 6.0; const UF: f64 = 0.026; const R: f64 = 1000.0; @@ -761,25 +820,25 @@ impl Samples { let mut system = System::new( ndim, |f: &mut Vector, x: f64, y: &Vector, _args: &mut NoArgs| { - let ue = C * f64::sin(D * x); - let f12 = BETA * (f64::exp((y[1] - y[2]) / UF) - 1.0); + let ue = A * f64::sin(OM * x); + let g12 = BETA * (f64::exp((y[1] - y[2]) / UF) - 1.0); f[0] = (y[0] - ue) / R; - f[1] = (2.0 * y[1] - UB) / S + GAMMA * f12; - f[2] = y[2] / S - f12; - f[3] = (y[3] - UB) / S + ALPHA * f12; + f[1] = (2.0 * y[1] - UB) / S + GAMMA * g12; + f[2] = y[2] / S - g12; + f[3] = (y[3] - UB) / S + ALPHA * g12; f[4] = y[4] / S; Ok(()) }, |jj: &mut CooMatrix, _x: f64, y: &Vector, m: f64, _args: &mut NoArgs| { - let g12 = BETA * f64::exp((y[1] - y[2]) / UF) / UF; + let h12 = BETA * f64::exp((y[1] - y[2]) / UF) / UF; jj.reset(); jj.put(0, 0, m * (1.0 / R)).unwrap(); - jj.put(1, 1, m * (2.0 / S + GAMMA * g12)).unwrap(); - jj.put(1, 2, m * (-GAMMA * g12)).unwrap(); - jj.put(2, 1, m * (-g12)).unwrap(); - jj.put(2, 2, m * (1.0 / S + g12)).unwrap(); - jj.put(3, 1, m * (ALPHA * g12)).unwrap(); - jj.put(3, 2, m * (-ALPHA * g12)).unwrap(); + jj.put(1, 1, m * (2.0 / S + GAMMA * h12)).unwrap(); + jj.put(1, 2, m * (-GAMMA * h12)).unwrap(); + jj.put(2, 1, m * (-h12)).unwrap(); + jj.put(2, 2, m * (1.0 / S + h12)).unwrap(); + jj.put(3, 1, m * (ALPHA * h12)).unwrap(); + jj.put(3, 2, m * (-ALPHA * h12)).unwrap(); jj.put(3, 3, m * (1.0 / S)).unwrap(); jj.put(4, 4, m * (1.0 / S)).unwrap(); Ok(())