(************** Content-type: application/mathematica ************** CreatedBy='Mathematica 5.2' Mathematica-Compatible Notebook This notebook can be used with any Mathematica-compatible application, such as Mathematica, MathReader or Publicon. The data for the notebook starts with the line containing stars above. To get the notebook into a Mathematica-compatible application, do one of the following: * Save the data starting with the line of stars above into a file with a name ending in .nb, then open the file inside the application; * Copy the data starting with the line of stars above to the clipboard, then use the Paste menu command inside the application. Data for notebooks contains only printable 7-bit ASCII and can be sent directly in email or through ftp in text mode. Newlines can be CR, LF or CRLF (Unix, Macintosh or MS-DOS style). NOTE: If you modify the data for this notebook not in a Mathematica- compatible application, you must delete the line below containing the word CacheID, otherwise Mathematica-compatible applications may try to use invalid cache data. For more information on notebooks and Mathematica-compatible applications, contact Wolfram Research: web: http://www.wolfram.com email: info@wolfram.com phone: +1-217-398-0700 (U.S.) Notebook reader applications are available free of charge from Wolfram Research. *******************************************************************) (*CacheID: 232*) (*NotebookFileLineBreakTest NotebookFileLineBreakTest*) (*NotebookOptionsPosition[ 26497, 759]*) (*NotebookOutlinePosition[ 27207, 784]*) (* CellTagsIndexPosition[ 27163, 780]*) (*WindowFrame->Normal*) Notebook[{ Cell[CellGroupData[{ Cell["Parametric Model Reduction Tutorial", "Title"], Cell[BoxData[ RowBox[{\(Evgenii\ B . \ Rudnyi\), ",", ButtonBox[\(\(\(http\)\(:\)\) // \(\(Evgenii . Rudnyi . Ru\)\(/\)\)\), ButtonData:>{ URL[ "http://Evgenii.Rudnyi.Ru/"], None}, ButtonStyle->"Hyperlink"]}]], "Text"], Cell[CellGroupData[{ Cell["Introduction", "Subtitle"], Cell[TextData[{ "This notebook allows you to repeat computations described in the paper\n\n\ E. B. Rudnyi, L. H. Feng, M. Salleras, S. Marco, J. G. Korvink. Error \ Indicator to Automatically Generate Dynamic Compact Parametric Thermal \ Models. THERMINIC 2005, 11th International Workshop on Thermal Investigations \ of ICs and Systems, 27 - 30 September 2005, Belgirate, Lake Maggiore, Italy, \ p. 139 - 145. ", ButtonBox["Preprint", ButtonData:>{ URL[ "http://www.imtek.uni-freiburg.de/simulation/mor4ansys/pdf/\ rudnyi05therminic.pdf"], None}, ButtonStyle->"Hyperlink"], " is at \ http://www.imtek.uni-freiburg.de/simulation/mor4ansys/pdf/rudnyi05therminic.\ pdf\n\nThe goal was to preserve three film coefficients for a simple thermal \ problem during model reduction process by means of multivariate expansion. \ The engineering problem is described in\n\nE. B. Rudnyi and J. G. Korvink. \ Boundary Condition Independent Thermal Model. In: Benner, P., Mehrmann, V., \ Sorensen, D. (eds) Dimension Reduction of Large-Scale Systems, Lecture Notes \ in Computational Science and Engineering (LNCSE). Springer-Verlag, \ Berlin/Heidelberg, Germany, v. 45, p. 345-348, 2005. ", ButtonBox["Preprint", ButtonData:>{ URL[ "http://www.imtek.uni-freiburg.de/simulation/benchmark/"], None}, ButtonStyle->"Hyperlink"], " is at http://www.imtek.uni-freiburg.de/simulation/benchmark/\n\nThe \ system matrices are available from ", ButtonBox["Oberwolfach Benchmark Collection", ButtonData:>{ URL[ "http://www.imtek.uni-freiburg.de/simulation/benchmark/"], None}, ButtonStyle->"Hyperlink"], " (for convenience they are included in the distribution along with this \ notebook)." }], "Text"], Cell["\<\ The original dynamic model after the discretization by the finite \ element method is as follows E dx/dt + (K + ht Kt + hs Ks + hb Kb) x = Bu The goal is to find a low-dimensional model in the similar form that \ approximates the dynamic properties of the original system in a wide range of \ film coefficients, ht, hs, hb. It is assumed that you browse the papers above before starting the tutorial. From a computational viewpoint, several Krylov subspaces are computed and \ then merged to form a projection basis. Several points in the parameter space \ of the transfer function are employed to choose an optimal dimension for each \ Krylov subspace. \ \>", "Text"] }, Open ]], Cell[CellGroupData[{ Cell["Setting up", "Subtitle"], Cell["\<\ Make the directory with system matrices current. You may need to \ change the next expression.\ \>", "Text"], Cell[BoxData[ \(SetDirectory["\<.\>"]\)], "Input"], Cell[TextData[{ "Load the functions. They are included for convenience but the latest \ versions as well as documentation are available from my ", ButtonBox["site", ButtonData:>{ URL[ "http://evgenii.rudnyi.ru/soft/Post4MOR.tar.gz"], None}, ButtonStyle->"Hyperlink"], " at http://evgenii.rudnyi.ru/soft/Post4MOR.tar.gz. The content of the \ archive is also available ", ButtonBox["online", ButtonData:>{ URL[ "http://evgenii.rudnyi.ru/soft/Post4MOR/"], None}, ButtonStyle->"Hyperlink"], " at http://evgenii.rudnyi.ru/soft/Post4MOR/. " }], "Text"], Cell[BoxData[ \(<< Post4MOR.m\)], "Input"], Cell[BoxData[ \(<< ModelReduction.m\)], "Input"], Cell["\<\ The first set of functions that is in the file Post4MOR.m is well \ documented. The goal of the second set was to add new functionality with the \ minimal development cost. As a result, my design was to reuse functions from \ Post4MOR.m as much as possible. Documentation for ModelReduction.m is in the \ beginning of the file.\ \>", "Text"], Cell[TextData[{ "Turn off some especially annoying ", StyleBox["Mathematica", FontSlant->"Italic"], " warning messages." }], "Text"], Cell[BoxData[{ \(\(Off[General::"\"];\)\), "\[IndentingNewLine]", \(\(Off[General::"\"];\)\)}], "Input"] }, Open ]], Cell[CellGroupData[{ Cell["Reading in and Preparing the Original System", "Subtitle"], Cell["\<\ The model in the benchmark was written as E dx/dt = (A - hbottom Abottom -htop Atop -hside Aside)x + Bu I convert it to the form used in the paper E dx/dt + (K + ht Kt + hs Ks + hb Kb) x = Bu K = -A Kt = Atop Ks = Aside Kb = Abottom Also, I leave only one output to reduce the number of plots. However, the \ procedure will work for many outputs as well.\ \>", "Text"], Cell[BoxData[{ \(\(matE\ = \ Import["\", "\"];\)\), "\[IndentingNewLine]", \(\(matK\ = \ \(-Import["\", "\"]\);\)\), "\ \[IndentingNewLine]", \(\(matB\ = \ Import["\", "\"];\)\), "\[IndentingNewLine]", \(\(matKt\ = \ Import["\", "\"];\)\), "\[IndentingNewLine]", \(\(matKs\ = \ Import["\", "\"];\)\), "\[IndentingNewLine]", \ \(\(matKb\ = \ Import["\", "\"];\)\), \ "\[IndentingNewLine]", \(\(matC\ = \ Take[Import["\", "\"], 1, Length[matK]];\)\)}], "Input"], Cell["\<\ Function MakeParametricSystem takes as an input a list of system \ matrices but we also have to define three functions to convert this list \ (mat) to three system matrices of DynamicSystem for a given list of \ parameters (par).\ \>", "Text"], Cell[BoxData[{ \(\(getM = Function[{mat, par}, Null];\)\), "\[IndentingNewLine]", \(\(getE = Function[{mat, par}, mat[\([1]\)]];\)\), "\[IndentingNewLine]", \(\(getK = Function[{mat, par}, mat[\([2]\)] + par[\([1]\)]*mat[\([3]\)] + par[\([2]\)]*mat[\([4]\)] + par[\([3]\)]*mat[\([5]\)]];\)\)}], "Input"], Cell[BoxData[ \(\(psys = MakeParametricSystem[{matE, matK, matKt, matKs, matKb}, matB, matC, {"\"}, {getM, getE, getK}];\)\)], "Input"], Cell["\<\ For better readability of the code below, I define several \ additional function to access system matrices from ParametricSystem.\ \>", \ "Text"], Cell[BoxData[{ \(MatrixE[sys_ParametricSystem] := Matrix[sys, 1]\), "\[IndentingNewLine]", \(MatrixK[sys_ParametricSystem] := Matrix[sys, 2]\), "\[IndentingNewLine]", \(MatrixKt[sys_ParametricSystem] := Matrix[sys, 3]\), "\[IndentingNewLine]", \(MatrixKs[sys_ParametricSystem] := Matrix[sys, 4]\), "\[IndentingNewLine]", \(MatrixKb[sys_ParametricSystem] := Matrix[sys, 5]\)}], "Input"] }, Open ]], Cell[CellGroupData[{ Cell["\<\ Computing Transfer Function for Several Points in the Parameter \ Space\ \>", "Subtitle"], Cell[TextData[{ "The question how to choose the points to check the local error is open. In \ the paper, the points along diagonals in the parameter space were used with \ an assumption that provided the approximation error is small for these points \ it is also small for internal points. \n\nWe compute\n\nH{s = 100, ht = 10, \ hb = 10, hs = 10}\nH{s = 100, ht = ", Cell[BoxData[ \(TraditionalForm\`10\^6\)]], ", hb = 10, hs = 10}\nH{s = 100, ht = ", Cell[BoxData[ \(TraditionalForm\`10\^6\)]], ", hb = ", Cell[BoxData[ \(TraditionalForm\`10\^6\)]], ", hs = 10}\nH{s = 100, ht = ", Cell[BoxData[ \(TraditionalForm\`10\^6\)]], ", hb = 10, hs = ", Cell[BoxData[ \(TraditionalForm\`10\^6\)]], "}" }], "Text"], Cell[TextData[{ "YSeries returns a matrix with the number of rows equal to the number of \ outputs and the number of columns equal to the number of simulation points \ (frequencies). As ", StyleBox["Mathematica", FontSlant->"Italic"], " stores matrices by rows and we need a column, I use Transpose. First \ takes the first row from the matrix." }], "Text"], Cell[BoxData[{ \(\(maxfreq = 100;\)\), "\[IndentingNewLine]", \(\(p1 = {10, 10, 10};\)\), "\[IndentingNewLine]", \(\(exact1\ = First[Transpose[ YSeries[\[IndentingNewLine]HarmonicSolution[{maxfreq}, MakeDynamicSystem[psys, p1]]\[IndentingNewLine]]]];\)\), "\[IndentingNewLine]", \(\(p2 = {1*^6, 10, 10};\)\), "\[IndentingNewLine]", \(\(exact2\ = \ First[Transpose[ YSeries[\[IndentingNewLine]HarmonicSolution[{maxfreq}, MakeDynamicSystem[psys, p2]]\[IndentingNewLine]]]];\)\), "\[IndentingNewLine]", \(\(p3 = {1*^6, 1*^6, 10};\)\), "\[IndentingNewLine]", \(\(exact3\ = \ First[Transpose[ YSeries[\[IndentingNewLine]HarmonicSolution[{maxfreq}, MakeDynamicSystem[psys, p3]]\[IndentingNewLine]]]];\)\), "\[IndentingNewLine]", \(\(p4 = {1*^6, 1*^6, 1*^6};\)\), "\[IndentingNewLine]", \(\(exact4\ = \ First[Transpose[ YSeries[\[IndentingNewLine]HarmonicSolution[{maxfreq}, MakeDynamicSystem[psys, p4]]\[IndentingNewLine]]]];\)\)}], "Input"], Cell["\<\ It could be good to take more points but in the general case, the \ computations with the original system are quite expensive.\ \>", "Text"] }, Open ]], Cell[CellGroupData[{ Cell["Preparing Common Things", "Subtitle"], Cell["\<\ An expansion point for the Laplace variable is zero, as this allows \ us to preserve the stationary solution. However, the expansion point for each \ parameter is 10. As a result, it is necessary to convert MatrixK to a new \ matrix.\ \>", "Text"], Cell["Transformation to hb0 = 10, hs0 = 10, hb0 = 10.", "Text"], Cell[BoxData[ \(matKnew\ = \ MatrixK[psys]\ + \ 10. *\((MatrixKt[psys] + MatrixKs[psys] + MatrixKb[psys])\)\)], "Input"], Cell["\<\ Note that the dimension of the original model is equal to \ 4257.\ \>", "Text"], Cell["\<\ Matrix is factorized only once. Later on this factor is used in \ order to quickly solve a system of equations matKnew x = b.\ \>", "Text"], Cell[BoxData[ \(\(invKnew = LinearSolve[matKnew, Method -> "\"];\)\)], "Input"], Cell["\<\ The starting vector will be the same for all Krylov \ subspaces.\ \>", "Text"], Cell[BoxData[ \(\(start = invKnew[MatrixB[psys] . {1}];\)\)], "Input"], Cell[TextData[{ "I multiply the input matrix by {1} in order to convert it from a matrix to \ a vector. In ", StyleBox["Mathematica", FontSlant->"Italic"], " a matrix consisting from a single column is different from a vector: a \ matrix 3x1 is {{1}, {2}, {3}} and a vector is {1, 2, 3}." }], "Text"], Cell["\<\ We will make a similar plots several times and below there are \ options to make it look nicer.\ \>", "Text"], Cell[BoxData[ \(\(plotoptions = {SymbolShape \[Rule] {PlotSymbol[Triangle, Filled \[Rule] False], PlotSymbol[Box, Filled \[Rule] False], PlotSymbol[Diamond, Filled \[Rule] False], PlotSymbol[Star, Filled \[Rule] False]}, SymbolStyle \[Rule] {RGBColor[1, 0, 0], RGBColor[0, 1, 0], RGBColor[0, 0, 1], RGBColor[0, 1, 1]}, FrameLabel \[Rule] {"\", "\<\>"}, PlotLegend \[Rule] {"\<10,10,10\>", "\<1e6,10,10\>", \ "\<1e6,1e6,10\>", "\<1e6,1e6,1e6\>"}, LegendBorder \[Rule] {}, LegendPosition \[Rule] {\(-0.7\), \(-0.4\)}, LegendSize \[Rule] {0.8, 0.28}, Evaluate[Imtek`Post4MOR`Private`defaultPlotOptions]};\)\)], "Input"] }, Open ]], Cell[CellGroupData[{ Cell["Generating Krylov Subspace Along the Laplace variable", "Subtitle"], Cell[TextData[{ "The function to compute the next vector is as follows \n\nnext = Knew^(-1) \ MatrixE prev \n\nand it is defined first. ArnoldiProcess uses this function \ and the starting vector defined above in order to compute the orthogonalized \ Krylov subspace \n\nKrylov{Knew^(-1) MatrixE v, Knew^(-1) f}\n\nThe \ projection matrix is returned in the transposed form as ", StyleBox["Mathematica", FontSlant->"Italic"], " stores matrices by rows. The maximum dimension 30 is taken based on \ previous experience." }], "Text"], Cell[BoxData[{ \(\(fun = Function[{prev}, invKnew[MatrixE[psys] . prev]];\)\), "\n", \(\(V = ArnoldiProcess[fun, start, 30];\)\), "\n", \(Dimensions[V]\)}], "Input"], Cell["Now we project a parametric system onto the basis V", "Text"], Cell[BoxData[ \(\(pred = ProjectSystem[psys, V];\)\)], "Input"], Cell["\<\ to make a reduced parametric system. Then we compute the relative \ error between the original system and the reduced system for the four \ different parameter sets defined above as a function of the system dimension. \ MakeDynamicSystem converts ParametricSystem to DynamicSystem for given set of \ parameters, FrequencyConvergence runs HarmonicSimulation for maxfreq for \ reduced systems with different dimensions from 1 to Length[MatrixK[red]]. \ Local error computes the relative error between results from reduced systems \ and the original system.\ \>", "Text"], Cell[BoxData[{ \(\(red1 = MakeDynamicSystem[pred, p1];\)\), "\[IndentingNewLine]", \(\(er1 = LocalError[FrequencyConvergence[maxfreq, red1], exact1];\)\), "\[IndentingNewLine]", \(\(red2 = MakeDynamicSystem[pred, p2];\)\), "\[IndentingNewLine]", \(\(er2 = LocalError[FrequencyConvergence[maxfreq, red2], exact2];\)\), "\[IndentingNewLine]", \(\(red3 = MakeDynamicSystem[pred, p3];\)\), "\[IndentingNewLine]", \(\(er3 = LocalError[FrequencyConvergence[maxfreq, red3], exact3];\)\), "\[IndentingNewLine]", \(\(red4 = MakeDynamicSystem[pred, p4];\)\), "\[IndentingNewLine]", \(\(er4 = LocalError[FrequencyConvergence[maxfreq, red4], exact4];\)\)}], "Input"], Cell[BoxData[ \(\(MultipleListPlot[{er1, er2, er3, er4}, Evaluate[plotoptions]];\)\)], "Input"], Cell["\<\ We see that the relative error for the set p1 decreases while the \ relative error for other sets (p2, p3, p4) stays high. This can be easily \ explained because we generate vectors for the Laplace variable and they do \ not approximate the system behavior well when film coefficients change \ considerably. In the paper, the local approximation error of 1% was considered to be good \ enough and the dimension 28 was chosen as the optimal for this \ subspace.\ \>", "Text"], Cell[BoxData[{ \(\(V\ = \ Take[V, \ 28];\)\), "\[IndentingNewLine]", \(Dimensions[V]\)}], "Input"] }, Open ]], Cell[CellGroupData[{ Cell["\<\ Generating Krylov Subspace Along the parameter ht (film coefficient \ at the top).\ \>", "Subtitle"], Cell["\<\ The function to compute the next vector now is as follows next = Knew^(-1) MatrixKt prev The maximum dimension 30 is taken rather arbitrary.\ \>", "Text"], Cell[BoxData[{ \(\(fun = Function[{prev}, invKnew[MatrixKt[psys] . prev]];\)\), "\n", \(\(W = ArnoldiProcess[fun, start, 30];\)\), "\n", \(Dimensions[W]\)}], "Input"], Cell["We merge subspaces V and W", "Text"], Cell[BoxData[{ \(\(V\ = \ Orthogonalize[V, W];\)\), "\[IndentingNewLine]", \(Dimensions[V]\)}], "Input"], Cell["\<\ One vector is deflated as the starting vector was the same for both \ subspaces.\ \>", "Text"], Cell["\<\ Now we project a parametric system onto the new basis V and compute \ relative errors again\ \>", "Text"], Cell[BoxData[ \(\(pred = ProjectSystem[psys, V];\)\)], "Input"], Cell[BoxData[{ \(\(red1 = MakeDynamicSystem[pred, p1];\)\), "\[IndentingNewLine]", \(\(er1 = LocalError[FrequencyConvergence[maxfreq, red1], exact1];\)\), "\[IndentingNewLine]", \(\(red2 = MakeDynamicSystem[pred, p2];\)\), "\[IndentingNewLine]", \(\(er2 = LocalError[FrequencyConvergence[maxfreq, red2], exact2];\)\), "\[IndentingNewLine]", \(\(red3 = MakeDynamicSystem[pred, p3];\)\), "\[IndentingNewLine]", \(\(er3 = LocalError[FrequencyConvergence[maxfreq, red3], exact3];\)\), "\[IndentingNewLine]", \(\(red4 = MakeDynamicSystem[pred, p4];\)\), "\[IndentingNewLine]", \(\(er4 = LocalError[FrequencyConvergence[maxfreq, red4], exact4];\)\)}], "Input"], Cell[BoxData[ \(\(MultipleListPlot[{er1, er2, er3, er4}, Evaluate[plotoptions]];\)\)], "Input"], Cell["\<\ We see that now the local error defined for the film coefficient ht \ goes down. Surprisingly, the local error for other two film coefficients goes \ down simultaneously. This can be explained that in our case study the film \ coefficient at the top plays the major role. Based on the convergence plot, we take the dimension of the reduced system of \ 41.\ \>", "Text"], Cell[BoxData[{ \(\(V\ = \ Take[V, \ 41];\)\), "\[IndentingNewLine]", \(Dimensions[V]\)}], "Input"] }, Open ]], Cell[CellGroupData[{ Cell["\<\ Generating Krylov Subspace Along the parameter hs (film coefficient \ at the side).\ \>", "Subtitle"], Cell["\<\ The function to compute the next vector now is as follows next = Knew^(-1) MatrixKs prev The maximum dimension 30 is taken rather arbitrary.\ \>", "Text"], Cell[BoxData[{ \(\(fun = Function[{prev}, invKnew[MatrixKs[psys] . prev]];\)\), "\n", \(\(W = ArnoldiProcess[fun, start, 30];\)\), "\n", \(Dimensions[W]\)}], "Input"], Cell["We merge subspaces V and W", "Text"], Cell[BoxData[{ \(\(V\ = \ Orthogonalize[V, W];\)\), "\[IndentingNewLine]", \(Dimensions[V]\)}], "Input"], Cell["\<\ One vector is deflated as the starting vector was the same for both \ subspaces.\ \>", "Text"], Cell["\<\ Now we project a parametric system onto the new basis V and compute \ relative errors again\ \>", "Text"], Cell[BoxData[ \(\(pred = ProjectSystem[psys, V];\)\)], "Input"], Cell[BoxData[{ \(\(red1 = MakeDynamicSystem[pred, p1];\)\), "\[IndentingNewLine]", \(\(er1 = LocalError[FrequencyConvergence[maxfreq, red1], exact1];\)\), "\[IndentingNewLine]", \(\(red2 = MakeDynamicSystem[pred, p2];\)\), "\[IndentingNewLine]", \(\(er2 = LocalError[FrequencyConvergence[maxfreq, red2], exact2];\)\), "\[IndentingNewLine]", \(\(red3 = MakeDynamicSystem[pred, p3];\)\), "\[IndentingNewLine]", \(\(er3 = LocalError[FrequencyConvergence[maxfreq, red3], exact3];\)\), "\[IndentingNewLine]", \(\(red4 = MakeDynamicSystem[pred, p4];\)\), "\[IndentingNewLine]", \(\(er4 = LocalError[FrequencyConvergence[maxfreq, red4], exact4];\)\)}], "Input"], Cell[BoxData[ \(\(MultipleListPlot[{er1, er2, er3, er4}, Evaluate[plotoptions]];\)\)], "Input"], Cell["\<\ We can see that the new subspace does not reduce the local error \ and we leave the dimension of the reduced system of 41.\ \>", "Text"], Cell[BoxData[{ \(\(V\ = \ Take[V, \ 41];\)\), "\[IndentingNewLine]", \(Dimensions[V]\)}], "Input"], Cell["\<\ The expansion along hb (the file coefficient at the bottom) \ produces the same result, that is, no further reduction in the local error. \ As a result, the basis of the dimension 41 is considered to be the optimal \ one. The optimal parametric reduced is system is\ \>", "Text"], Cell[BoxData[ \(\(pred = ProjectSystem[psys, V];\)\)], "Input"] }, Open ]], Cell[CellGroupData[{ Cell["Comparing the Reduced and Original systems", "Subtitle"], Cell["\<\ The procedure above is empirical by its nature. There were two main \ assumptions made: 1) One can ignore mixed moments. 2) The relative error at a few border points in the parameter space is enough \ to choose the dimension of the reduced model. In our experience, this was working okay for the problem in question but it \ is unclear whether this will work for other systems. Unfortunately, \ mathematical proofs to estimate errors in the case of parametric model \ reduction are missing. Below, an empirical check of the quality of the reduced system is \ performed.\ \>", "Text"], Cell["\<\ Function that compares the static, harmonic response and transient \ solution of the original and reduced systems. We compute relative error in \ per cent as follows Norm[full-red]/Norm[full]*100 for each output and then \ the maximum values is taken to the final table.\ \>", "Text"], Cell[BoxData[{ \(error[r2_SimulationResult, r1_SimulationResult] := MapThread[ Norm[#1 - #2]/Norm[#2] &, {YSeries[r2], YSeries[r1]}]\[IndentingNewLine]\), "\[IndentingNewLine]", \(compareSystems[psys_ParametricSystem, pred_ParametricSystem, par_List] := Module[{sys, red, r1, r2, relS, relT, relH, time = Table[10^i, {i, \(-5\), 4, 0.1}], freq = Table[10^i, {i, \(-4\), 3, 0.5}]}, \[IndentingNewLine]Print[ par]; \[IndentingNewLine]sys = MakeDynamicSystem[psys, par]; \[IndentingNewLine]red = MakeDynamicSystem[pred, par]; \[IndentingNewLine]r1 = StationarySolution[red]; \[IndentingNewLine]r2 = StationarySolution[sys]; \[IndentingNewLine]relS = Abs[r1 - r2]/ r2*100; \[IndentingNewLine]Print["\", relS]; \[IndentingNewLine]r1 = AnsysTransientSolution[time, red]; \[IndentingNewLine]r2 = AnsysTransientSolution[time, sys]; \[IndentingNewLine]PlotResult[{r2, r1}, PlotStyle \[Rule] {RGBColor[1, 0, 0], RGBColor[0, 1, 0]}, FunctionX \[Rule] Log10, CommonTitle -> "\"]; \ \[IndentingNewLine]PlotResult[ Difference[r2, r1, ErrorFunction \[Rule] \((Log10[Abs[#1 - #2]/Abs[#1]] &)\)], PlotStyle \[Rule] {RGBColor[0, 1, 0]}, FunctionX \[Rule] Log10, CommonTitle -> "\"]; \ \[IndentingNewLine]relT = error[r2, r1]*100; \[IndentingNewLine]Print["\", relT]; \[IndentingNewLine]r1 = HarmonicSolution[freq, red]; r2 = HarmonicSolution[freq, sys]; \[IndentingNewLine]PlotResult[{r2, r1}, FunctionY \[Rule] \((Log10[Abs[#]] &)\), FunctionX \[Rule] Log10, PlotStyle \[Rule] {RGBColor[1, 0, 0], RGBColor[0, 1, 0]}, CommonTitle -> "\"]; PlotResult[ Difference[r2, r1, ErrorFunction \[Rule] \((Log10[Abs[#1 - #2]/Abs[#1]] &)\)], FunctionX \[Rule] Log10, PlotStyle \[Rule] {RGBColor[0, 1, 0]}, CommonTitle -> "\"]; relH = error[r2, r1]*100; \[IndentingNewLine]Print["\", relH]; \[IndentingNewLine]Flatten[{par, Max[relS], Max[relT], Max[relH]}]]\)}], "Input"], Cell["\<\ Below there is a check the function above for some parameter \ values. It takes some time as it makes simulation for the original system as \ well. First plot contains two curves: a red line for the original system and \ a green line for the reduced system. As the difference is small, one can see \ red just a little bit. The second plot gives the relative difference between \ the two curves.\ \>", "Text"], Cell[BoxData[ \(compareSystems[psys, pred, {1, 100, 1000}]\)], "Input"], Cell["\<\ Now we make a loop over different combinations of parameters. \ Computational time is equal to that above multiplied by \ Length[krange]^3.\ \>", "Text"], Cell[BoxData[ \(\(krange\ = \ {1, \ 10000};\)\)], "Input"], Cell[BoxData[ \(\(table\ = Outer[compareSystems[psys, pred, {#1, #2, #3}] &, \ krange, \ krange, \ krange];\)\)], "Input"], Cell[BoxData[ \(TableForm[Flatten[table, 2]]\)], "Input"], Cell["\<\ We can see that the low-dimensional system of dimension 41 does \ approximate the behavior of the original system of dimension 4257 quite well.\ \ \>", "Text"] }, Open ]], Cell[CellGroupData[{ Cell[TextData[{ "Bringing ", StyleBox["Mathematica", FontSlant->"Italic"], " Properties into the Default State" }], "Subtitle"], Cell[BoxData[{ \(\(On[General::"\"];\)\), "\[IndentingNewLine]", \(\(On[General::"\"];\)\)}], "Input"] }, Open ]] }, Open ]] }, FrontEndVersion->"5.2 for Macintosh", ScreenRectangle->{{0, 1024}, {0, 702}}, WindowSize->{739, 585}, WindowMargins->{{Automatic, 5}, {Automatic, 0}}, PrintingCopies->1, PrintingPageRange->{1, Automatic}, ShowSelection->True ] (******************************************************************* Cached data follows. If you edit this Notebook file directly, not using Mathematica, you must remove the line containing CacheID at the top of the file. The cache data will then be recreated when you save this file from within Mathematica. *******************************************************************) (*CellTagsOutline CellTagsIndex->{} *) (*CellTagsIndex CellTagsIndex->{} *) (*NotebookFileOutline Notebook[{ Cell[CellGroupData[{ Cell[1776, 53, 52, 0, 168, "Title"], Cell[1831, 55, 255, 5, 35, "Text"], Cell[CellGroupData[{ Cell[2111, 64, 32, 0, 59, "Subtitle"], Cell[2146, 66, 1742, 34, 359, "Text"], Cell[3891, 102, 682, 16, 245, "Text"] }, Open ]], Cell[CellGroupData[{ Cell[4610, 123, 30, 0, 59, "Subtitle"], Cell[4643, 125, 118, 3, 36, "Text"], Cell[4764, 130, 54, 1, 33, "Input"], Cell[4821, 133, 584, 14, 74, "Text"], Cell[5408, 149, 46, 1, 33, "Input"], Cell[5457, 152, 52, 1, 33, "Input"], Cell[5512, 155, 351, 6, 93, "Text"], Cell[5866, 163, 142, 5, 36, "Text"], Cell[6011, 170, 130, 2, 52, "Input"] }, Open ]], Cell[CellGroupData[{ Cell[6178, 177, 64, 0, 59, "Subtitle"], Cell[6245, 179, 382, 15, 283, "Text"], Cell[6630, 196, 718, 17, 147, "Input"], Cell[7351, 215, 253, 5, 74, "Text"], Cell[7607, 222, 372, 8, 109, "Input"], Cell[7982, 232, 168, 3, 52, "Input"], Cell[8153, 237, 155, 4, 55, "Text"], Cell[8311, 243, 437, 9, 109, "Input"] }, Open ]], Cell[CellGroupData[{ Cell[8785, 257, 99, 3, 91, "Subtitle"], Cell[8887, 262, 759, 21, 210, "Text"], Cell[9649, 285, 367, 8, 74, "Text"], Cell[10019, 295, 1211, 25, 337, "Input"], Cell[11233, 322, 150, 3, 55, "Text"] }, Open ]], Cell[CellGroupData[{ Cell[11420, 330, 43, 0, 59, "Subtitle"], Cell[11466, 332, 257, 5, 74, "Text"], Cell[11726, 339, 63, 0, 36, "Text"], Cell[11792, 341, 157, 4, 71, "Input"], Cell[11952, 347, 89, 3, 36, "Text"], Cell[12044, 352, 149, 3, 55, "Text"], Cell[12196, 357, 103, 2, 33, "Input"], Cell[12302, 361, 88, 3, 36, "Text"], Cell[12393, 366, 74, 1, 33, "Input"], Cell[12470, 369, 310, 7, 74, "Text"], Cell[12783, 378, 119, 3, 36, "Text"], Cell[12905, 383, 766, 12, 242, "Input"] }, Open ]], Cell[CellGroupData[{ Cell[13708, 400, 73, 0, 91, "Subtitle"], Cell[13784, 402, 541, 10, 226, "Text"], Cell[14328, 414, 179, 3, 71, "Input"], Cell[14510, 419, 67, 0, 36, "Text"], Cell[14580, 421, 67, 1, 33, "Input"], Cell[14650, 424, 578, 9, 131, "Text"], Cell[15231, 435, 773, 16, 166, "Input"], Cell[16007, 453, 108, 2, 33, "Input"], Cell[16118, 457, 486, 10, 150, "Text"], Cell[16607, 469, 108, 2, 52, "Input"] }, Open ]], Cell[CellGroupData[{ Cell[16752, 476, 110, 3, 91, "Subtitle"], Cell[16865, 481, 167, 6, 112, "Text"], Cell[17035, 489, 180, 3, 71, "Input"], Cell[17218, 494, 42, 0, 36, "Text"], Cell[17263, 496, 114, 2, 52, "Input"], Cell[17380, 500, 104, 3, 36, "Text"], Cell[17487, 505, 115, 3, 36, "Text"], Cell[17605, 510, 67, 1, 33, "Input"], Cell[17675, 513, 773, 16, 166, "Input"], Cell[18451, 531, 108, 2, 33, "Input"], Cell[18562, 535, 381, 8, 112, "Text"], Cell[18946, 545, 108, 2, 52, "Input"] }, Open ]], Cell[CellGroupData[{ Cell[19091, 552, 111, 3, 91, "Subtitle"], Cell[19205, 557, 167, 6, 112, "Text"], Cell[19375, 565, 180, 3, 71, "Input"], Cell[19558, 570, 42, 0, 36, "Text"], Cell[19603, 572, 114, 2, 52, "Input"], Cell[19720, 576, 104, 3, 36, "Text"], Cell[19827, 581, 115, 3, 36, "Text"], Cell[19945, 586, 67, 1, 33, "Input"], Cell[20015, 589, 773, 16, 166, "Input"], Cell[20791, 607, 108, 2, 33, "Input"], Cell[20902, 611, 146, 3, 55, "Text"], Cell[21051, 616, 108, 2, 52, "Input"], Cell[21162, 620, 289, 5, 74, "Text"], Cell[21454, 627, 67, 1, 33, "Input"] }, Open ]], Cell[CellGroupData[{ Cell[21558, 633, 62, 0, 59, "Subtitle"], Cell[21623, 635, 596, 15, 226, "Text"], Cell[22222, 652, 294, 5, 74, "Text"], Cell[22519, 659, 2528, 47, 755, "Input"], Cell[25050, 708, 418, 7, 93, "Text"], Cell[25471, 717, 75, 1, 33, "Input"], Cell[25549, 720, 163, 4, 55, "Text"], Cell[25715, 726, 63, 1, 33, "Input"], Cell[25781, 729, 147, 3, 52, "Input"], Cell[25931, 734, 61, 1, 33, "Input"], Cell[25995, 737, 169, 4, 55, "Text"] }, Open ]], Cell[CellGroupData[{ Cell[26201, 746, 137, 5, 91, "Subtitle"], Cell[26341, 753, 128, 2, 52, "Input"] }, Open ]] }, Open ]] } ] *) (******************************************************************* End of Mathematica Notebook file. *******************************************************************)