Re: [Maxima-discuss] Construct sequence of identifiers

[Robert D. <rob...@gm...>]

Wolfgang, here is an implementation of the finite difference function which
makes use of subscripted variables, as I was suggesting.

I have arranged the main function so that it only constructs the list of
equations and does not solve them -- I think it's important to review the
equations before solving. Also, I'm using centered differences since that
way the indexing is easier for me to keep straight.

By the way, the functions x, yp, and ypp are not really local functions;
all functions defined by ":=" are global. But making them actually local
requires introducing lambda expressions, which is perhaps beyond the scope
of this.

finite_difference_eqs (x_0, x_n, y_0, y_n, n) :=
    block ([h: (x_n - x_0) / n],
           x(k) := x_0 + k*h,
           yp(k) := (y[k + 1] - y[k - 1]) / (2*h),
           ypp(k) := (y[k + 1] - 2*y[k] + y[k - 1])/h^2,
           append ([y[0] = y_0], makelist (BVP(i), i, 1, n - 1), [y[n] =
y_n]));

BVP(k) := ypp(k) + (2/10)*yp(k) + 4*y[k] = 3*x(k) - 1;

Here's the output I get for the problem which you stated.

(%i5) n:4
(%o5)                           4
(%i6) eqs:finite_difference_eqs(0,1,1/10,7/10,n)
                                      2 (y  - y )
            1                             2    0             1
(%o6) [y  = --, 16 (y  - 2 y  + y ) + ----------- + 4 y  = - -,
        0   10       2      1    0         5           1     4
                      2 (y  - y )
                          3    1           1
16 (y  - 2 y  + y ) + ----------- + 4 y  = -,
     3      2    1         5           2   2
                      2 (y  - y )
                          4    2           5       7
16 (y  - 2 y  + y ) + ----------- + 4 y  = -, y  = --]
     4      3    2         5           3   4   4   10

Okay, now we have the list of equations. Now solve the equations.

(%i7) yy:makelist(y[i],i,0,n)
(%o7)                 [y , y , y , y , y ]
                        0   1   2   3   4
(%i8) linsolve(eqs,yy)
            1        543411        22751       878931        7
(%o8) [y  = --, y  = -------, y  = -----, y  = -------, y  = --]
        0   10   1   1191400   2   34040   3   1191400   4   10

With the list of equations in hand, we can inspect it in different ways.
Here I'll look at the augmented matrix of coefficients.

(%i9) augcoefmatrix(eqs,yy)
               [                             1  ]
               [ 1    0     0     0    0   - -- ]
               [                             10 ]
               [                                ]
               [ 78         82              1   ]
               [ --  - 28   --    0    0    -   ]
               [ 5          5               4   ]
               [                                ]
               [      78          82         1  ]
(%o9)          [ 0    --   - 28   --   0   - -  ]
               [      5           5          2  ]
               [                                ]
               [            78         82    5  ]
               [ 0    0     --   - 28  --  - -  ]
               [            5          5     4  ]
               [                                ]
               [                             7  ]
               [ 0    0     0     0    1   - -- ]
               [                             10 ]


Hope this helps,

Robert