Re: [Maxima-discuss] Construct sequence of identifiers

[Dr. W. L. <dr....@gm...>]

Dear Robert,

This is very good news.
I was thinking about yours and Richard suggestions ..
And now I can learn much from your code :) Wonderful.

Thank you for your kind lesson, I like it.

Best regards with sympathetic feelings for Maxima,
Wolfgang

> Am 21.07.2023 um 21:17 schrieb Robert Dodier <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
>