Quit[]

\(\(\( (*Solve\ non - 
        homogenouse\ first\ order\ ODE\ with\ constant\ coeffs\ using\ the\ FEM\ method*) \)\(
\)\( (*by\ Nasser\ Abbasi . \ 
          Sept\ 2006*) \)\(
\)\( (*\ 
      a\ dy/dt\  + \ b\ y \((t)\)\  = \ f \((t)\), \ 
      with\ initial\ condition\ y \((t0)\) = 
        y0\ *) \)\(
\)\(
\)\(<< 
      LinearAlgebra`MatrixManipulation`

    << Graphics`

    \(SetDirectory["E:/nabbasi/data/nabbasi_web_Page/my_courses/UCI_COURSES/CREDIT_COURSES/spring_2006/spring_MAE_207/other_class_documents/Nasser_Abbasi/small_report"];\)

 (*********************************************************) 

    getStiffnessMatrix[a_, b_, h_, nNodes_] := 
      Module[{e1, e2, e3, m = nNodes - 1, k}, 
k = 
          Table[0, {m}, {m}]; 
e1 = \(2.  b\ h\)\/3; \ 
        e2 = \(3.  a + b\ h\)\/6; e3 = \(\(-3. \) a + b\ h\)\/6; 
        e4 = \(b\ h - 3.  a\)\/6; 
        e5 = \(3.  a + 2  b\ h\)\/6; 
k[\([1, 1]\)] = e1; \ 
        k[\([1, 2]\)] = e2; 
k[\([m, m - 1]\)] = e4; \ 
        k[\([m, m]\)] = e5; 
For[i = 2, 
          i ≤ 
            m - 1, \(i++\), 
k[\([i, i - 1]\)] = 
            e3; 
k[\([i, i]\)] = 
            e1; 
k[\([i, i + 1]\)] = 
            e2;
]; 
k
]

 (*********************************************************) 

    getLoadVectorOld[a_, b_, h_, f_, u0_, nNodes_, length_, xc_] := 
      Module[{m = nNodes - 1, load, 
          nElements = nNodes - 1}, 
load = 
          Table[0, {i, 1, 
              m}]; 

load[\([1]\)] = \(1. \/h\) \((∫\_\(xc[\([1]\)]\)\%\(xc[\([2]\)]\)\((\ \((\ 
                            x - xc[\([1]\)]\ )\) 
                          f[x]\ )\) \[DifferentialD]x + ∫\_\(xc[\([2]\)]\)\%\(xc[\([3]\)]\)\((\ \((\ xc[\([3]\)] - x\ )\) 
                          f[
                            x]\ )\) \[DifferentialD]x)\) - \(\((\(-3\) a + b\ h)\) u0\)\/6; 
load[\([m]\)] = \(1. \/h\) \(∫\_\(xc[\([nNodes - 1]\)]\)\%\(xc[\([nNodes]\)]\)\((\((x - 
                        xc[\([nNodes - 1]\)])\) 
                    f[x])\) \[DifferentialD]x\); 

For[i = 2, 
          i ≤ 
            m - 1, \(i++\), 
\(load[\([i]\)] = \(1. \/h\) \((∫\_\(xc[\([i]\)]\)\%\(xc[\([i + 1]\)]\)\((\((x - xc[\([i]\)])\) 
                            f[x])\) \[DifferentialD]x + ∫\_\(xc[\([i + 1]\)]\)\%\(xc[\([i + 2]\)]\)\((\((xc[\([i + 2]\)] - 
                                x)\) f[
                              x])\) \[DifferentialD]x)\);\)
]; 
load
]

 (*********************************************************) 

    getLoadVector[a_, b_, h_, f_, u0_, nNodes_, length_, xc_] := 
      Module[{m = nNodes - 1, load, 
          nElements = nNodes - 1}, 
load = 
          Table[0, {i, 1, 
              m}]; 

load[\([1]\)] = \(1\/h\) NIntegrate[\((x - xc[\([1]\)])\) f[x], {x, xc[\([1]\)], 
                  xc[\([2]\)]}] + \(1\/h\) 
              NIntegrate[\((xc[\([3]\)] - x)\) f[x], {x, xc[\([2]\)], 
                  xc[\([3]\)]}] - \(\((\(-3\) a + b\ h)\) u0\)\/6; 
load[\([m]\)] = \(1\/h\) 
            NIntegrate[\((x - xc[\([nNodes - 1]\)])\) f[x], {x, 
                xc[\([nNodes - 1]\)], 
                xc[\([nNodes]\)]}]; 

For[
          i = 2, i ≤ 
            m - 1, \(i++\), 
load[\([i]\)] = \(1\/h\) 
                NIntegrate[\((x - xc[\([i]\)])\) f[x], {x, xc[\([i]\)], 
                    xc[\([i + 1]\)]}] + \(1\/h\) 
                NIntegrate[\((xc[\([i + 2]\)] - x)\) f[x], {x, 
                    xc[\([i + 1]\)], 
                    xc[\([i + 
                          2]\)]}]
]; 
load
]

 (*********************************************************) 

    getRmsError[h_, nNodes_, xc_, x0_, analyticalSolution_, length_, q_, 
        theTrialSolution_, x_] := 
      Module[{delX, nPoints, rmserror, sum, t, v, 
          at}, 
delX = length\/1000. ; \ 
        nPoints = 1001; 

sum = 
          0. ; 
For[i = 0, 
          i ≤ nPoints - 1, \(i++\), 
at = 
            x0 + i*delX; 
t = 
            analyticalSolution /. x -> at; 
v = 
            theTrialSolution /. x -> at; 
sum = 
            sum + \((t - v)\)\^2;
]; 
rmserror = \@\(\(\ \)\(sum\/nPoints\)\); 
rmserror = 
          rmserror[\([1]\)]
]

 (*********************************************************) 

    getSolution[a_, b_, h_, nNodes_, x0_, f_, u0_, length_, xc_] := 
      Module[{q, k, load, sol, z}, 

q = 
          Table[0, {i, 1, nNodes}]; 

z = 
          Timing[getStiffnessMatrix[a, b, h, nNodes]]; 
k = 
          z[\([2]\)]; 
Print["Timing> getStiffnessMatrix=", z[\([1]\)]]; 

z = 
          Timing[getLoadVector[a, b, h, f, u0, nNodes, length, 
              xc]]; 
load = 
          z[\([2]\)]; 
Print["Timing> getloadVector=", 
          z[\([1]\)]]; 

sol = 
          LinearSolve[k, 
            load]; 
q[\([\ Range[2, nNodes]\ ]\)] = 
          sol; 
q[\([1]\)] = 
          u0; 
q
]

 (*********************************************************) 

    fem[a_, b_, nNodes_, x0_, f_, u0_, length_, rhsString_, miny_, maxy_] := 
      Module[{q, eq, xc, pts, p, p1, p2, h, nPoints, delX, rmserror, 
          theTrialSolution, x}, 
        xc = Table[i  length\/\(\(nNodes - 1\)\(
\)
                \), {i, 0, nNodes - 1}]; 
xc = 
          xc + x0; 
h = 
          xc[\([2]\)] - xc[\([1]\)]; 
q = 
          getSolution[a, b, h, nNodes, x0, f, u0, length, 
            xc]; 
eq = 
          a\ \(u'\)[x] + b\ u[x] == f[x]; 
sol1 = 
          Chop[DSolve[{\ \ eq, u[x0] == u0}, \ u[x], 
              x]]; 
analyticalSolution = 
          u[x] /. sol1; 

theTrialSolution = 
          trial[x, h, q, nNodes, 
            xc]; 

z = 
          Timing[getRmsError[h, nNodes, xc, x0, analyticalSolution, length, 
              q, theTrialSolution, x]]; 
rmserror = 
          z[\([2]\)]; 
Print["Timing> getRmsError=", 
          z[\([1]\)]]; 

title = 
          If[a ≠ 1, ToString[a], ""] <> \ "y'[x] +" <> 
            If[b ≠ 1, ToString[b], ""] <> " y[x]=" <> 
            rhsString <> " ,y(0)=" <> ToString[u0] <> "
RMS ERROR=" <> 
            ToString[N[rmserror]] <> "
number of nodes=" <> 
            ToString[
              nNodes]\ ; 

 (*\(title = 
              Append[{title}, {"
Analytical solution= "\ , 
                  analyticalSolution}];\)*) 
 (*\(title = 
              Flatten[title];\)*) 

z = 
          Timing[
Plot[{theTrialSolution, 
                analyticalSolution}, {x, 0, length}, 
              PlotRange -> {{0, length + 0.01\ length}, {miny, maxy}}, 
              AxesOrigin -> {0, 0}, 
              PlotStyle -> {Hue[0.9], Hue[ .6]}, 
              PlotLabel -> title, 
              TextStyle -> {FontSize -> 16}, 
              DisplayFunction -> 
                Identity
]
]; 
p1 = 
          z[\([2]\)]; 
Print["Timing> plot first=", 
          z[\([1]\)]]; 

pts = 
          Table[{xc[\([i]\)], theTrialSolution /. x -> xc[\([i]\)]}, {i, 
              1, nNodes}]; 
z = 
          Timing[ListPlot[pts, PlotStyle -> PointSize[0.02], 
              PlotRange -> {{0, length + 0.01\ length}, {miny, maxy}}, 
              AxesOrigin -> {0, 0}, 
              DisplayFunction -> 
                Identity]
]; 
p2 = 
          z[\([2]\)]; 
Print["Timing> plot second=", 
          z[\([1]\)]]; 

 (*\(showShapeFunctions[];\)*) 
p = 
          Show[{p1, p2}, 
            DisplayFunction -> 
              Identity]
]

 (*********************************************************) 

    trial[x_, h_, q_, nNodes_, xc_] := 
      Module[{}, 
Sum[
          q[\([i]\)]\ phi[x, i, h, nNodes, xc], {i, 
            nNodes}]
]

 (*********************************************************) 

    N1[i_, h_, x_, xc_] := 
      Module[{}, \(xc[\([i + 1]\)] - x\)\/h]

 (*********************************************************) 

    N2[i_, h_, x_, xc_] := 
      Module[{}, \(x - xc[\([i - 1]\)]\)\/h]

 (*********************************************************) 

    phi[x_, j_, h_, nNodes_, xc_] := 
      Module[{r}, 

r = 
          Table[Which[i == 1, 
\ 
              N1[i, h, x, 
                  xc] \((UnitStep[x - xc[\([i]\)]] - 
                    UnitStep[
                      x - xc[\([i + 1]\)]])\), 
i == 
                nNodes, 
              N2[i, h, x, 
                  xc] \((UnitStep[
                    x - xc[\([i - 1]\)]])\), 
True, 
              N2[i, h, x, 
                    xc] \((UnitStep[x - xc[\([i - 1]\)]] - 
                      UnitStep[x - xc[\([i]\)]])\) + 
                N1[i, h, x, 
                    xc] \((UnitStep[x - xc[\([i]\)]] - 
                      UnitStep[
                        x - xc[\([i + 
                                1]\)]])\)
], 
{i, 1, 
              nNodes}
]; 

r[\([j]\)]
]

    \)\)\)

\(\(\( (*********************************************************) \)\(
\)\( (*The\ MAIN\ starts\ here\ *) \)\(
\)\(
\)\( (*\ 
      solves\ \ \ \ a\ dy/dt\  + \ b\ y \((t)\)\  = \ 
        rhs\ \ \ with\ initial\ condition\ y0\ *) \)\(
\)\(
\)\(\(fileit[fileName_, p_] := 
      Module[{}, 
\(Export[
            fileName <> ".gif", \ 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ p, \ 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ ConversionOptions\  -> \ {"AnimationDisplayTime"\  -> \  .5, \ "Loop"\  -> \ 
                  False, "GlobalColorReduction" -> 
                  True}, \ 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ ImageSize\  -> \ {640, \ 
                Automatic}
\ \ \ \ \ \ \ \ \ \ \ \ ];\)
]\)\(
\)
    \)\)\)

\(\(\( (*\ 
      solves\ \ \ \ \ dy/dt\  + \ 2\ y \((t)\)\  = \ 
        Sin[x]\ \ \ with\ initial\ condition\ y0\ *) \)\(
\)\(length = 15. ; \ maxNodes = 100; \ a = 1; \ 
    b = 2; \ \ \ p = {}; \ x0 = 0; \ u0 = 1; \ miny = \(- .7\); 
    maxy = 1.1;

    f[x_] := Module[{}, Sin[x]]

    \(Do[Print["i=", i]; \ 
        p = Append[p, 
            fem[a, b, i, x0, f, u0, length, "sin(x)", miny, 
              maxy]];\ 
, {i, 3, 
          maxNodes}];\)
 (*Do[
        Show[p[\([i]\)], DisplayFunction -> $DisplayFunction], {i, 1, 
          Length[p]}]*) 

    \(fileit["eq1", p];\)\)\)\)

\(\(\(
\)\( (*\ 
      solves\ \ \ \ \ dy/dt\  + \  .5\ y \((t)\)\  = \ \(cos \((x)\) 
            Sin[x]\ \ \ with\ initial\ condition\ y0 = \(-1\)\)\ *) \)\(
\)\(length = 10. ; \ maxNodes = 50; \ a = 1; \ 
    b =  .2; \ \ p = {}; \ x0 = 0; \ u0 = \(-1\); \ miny = \(-1.1\); 
    maxy =  .5;

    f[x_] := Module[{}, Sin[x]\ Cos[x]]

    \(Do[Print["i=", i]; \ 
        p = Append[p, 
            fem[a, b, i, x0, f, u0, length, "cos(x) sin(x)", miny, 
              maxy]];\ 
, {i, 3, 
          maxNodes}];\)
 (*Do[
        Show[p[\([i]\)], DisplayFunction -> $DisplayFunction], {i, 1, 
          Length[p]}]*) 

    \(fileit["eq2", p];\)\)\)\)

\(\(\(
\)\( (*\ 
      solves\ \ \ \ \ 5\ dy/dt\  + \ 
          20\ y \((t)\)\  = \ \(x^3\ Sin[
              x]\ \ \ with\ initial\ condition\ y0 =  .6\)\ *) \)\(
\)\(length = 4. ; \ maxNodes = 100; \ a = 5; \ 
    b = 20; \ \ p = {}; \ x0 = 0; \ u0 =  .6; \ miny = \(- .1\); 
    maxy =  .6;

    f[x_] := Module[{}, x^3\ Sin[x]]

    \(Do[Print["i=", i]; \ 
        p = Append[p, 
            fem[a, b, i, x0, f, u0, length, "x^3 * sin(x)", miny, 
              maxy]];\ 
, {i, 3, 
          maxNodes}];\)
 (*Do[
        Show[p[\([i]\)], DisplayFunction -> $DisplayFunction], {i, 1, 
          Length[p]}]*) 

    \(fileit["eq3", p];\)\)\)\)

\(\(\(
\)\( (*\ 
      solves\ \ \ \ \ 2\ dy/dt\  + \ 
          10\ y \((t)\)\  = \ \(x^4\ Cos[x]\ Sin[
              x]\ \ \ with\ initial\ condition\ y0 = 
          6\)\ *) \)\(
\)\(length = 6; \ maxNodes = 100; \ 
    a = 2; \ b = 10; \ \ p = {}; \ x0 = 0; \ u0 = 6; \ miny = \(-10\); 
    maxy = 14;

    f[x_] := Module[{}, x^4\ Cos[x]\ Sin[x]]

    \(Do[Print["i=", i]; \ 
        p = Append[p, 
            fem[a, b, i, x0, f, u0, length, "x^3 * sin(x)", miny, 
              maxy]];\ 
, {i, 3, 
          maxNodes}];\)
 (*Do[
        Show[p[\([i]\)], DisplayFunction -> $DisplayFunction], {i, 1, 
          Length[p]}]*) 

    \(fileit["eq4", p];\)

    \)\)\)