First, let us set display2d to false. This is suitable for writing code into e-mails. On the other
hand the matrices and a lot of other things aren't printed pretty then. It will look much nicer
when you omit the line display2d : false on your computer.
(%i1) display2d : false$
Let's create now a simple linear system with a simple solution:
(%i2) 2*x+3*y+4*z, [x:1,y:2,z:3];
(%o2) 20
I asked, what will 2*x+3*y+4*z be if I locally set x to 1, etc. (In top level x is still x.)
(%i3) 3*x-y-3*z, [x:1,y:2,z:3];
(%o3) -8
(%i4) -x-y+4*z, [x:1,y:2,z:3];
(%o4) 9
OK, here is the linear system:
(%i5) sys: [2*x+3*y+4*z-20=0, 3*x-y-3*z+8=0, -x-y+4*z-9=0]$
Let A be the augmented coefficient matrix:
(%i6) A: augcoefmatrix (sys, [x,y,z]);
(%o6) matrix([2,3,4,-20],[3,-1,-3,8],[-1,-1,4,-9])
What you want I guess is done by the function echelon:
(%i7) echelon(A);
(%o7) matrix([1,3/2,2,-10],[0,1,18/11,-76/11],[0,0,1,-3])
(It really will look much nicer if you omit display2d : false.)
Before we start to implement your pseudocode, you should know some things about
matrices in Maxima:
(%i8) col(A,1);
(%o8) matrix([2],[3],[-1])
(%i9) row(A,1);
(%o9) matrix([2,3,4,-20])
col and row still return matrices, while args returns a list of lists of row entries and A[i]
returns row i as a list.
(%i10) a: args(A);
(%o10) [[2,3,4,-20],[3,-1,-3,8],[-1,-1,4,-9]]
(%i11) A[1];
(%o11) [2,3,4,-20]
We can calculate the number m of rows and number n of cols (coefficients):
(%i12) m: length(a);
(%o12) 3
(%i13) n: length(A[1])-1;
(%o13) 3
Now let's define a function which implements your pseudocode:
gaussian_elimination (A) := block( [i, j, m, n, maxi, listarith:true],
...
A )$
The function gaussian_elimination takes A as argument and returns A which will be newly
calculated by your code. i, j, m, n, maxi are defined to be local variables inside of block.
Also I set listarith to true (which is the default, but I need it --> see below) locally.
This is the skeleton. Now let's fill it with your code:
gaussian_elimination (A) := block( [i, j, m, n, maxi, listarith:true],
m:length(args(A)),
n:length(A[1])-1,
i:1,
j:1,
while (i<=m and j<=n) do (
maxi:i,
for k:i+1 thru m do
if abs(A[k,j]) > abs(A[maxi,j]) then
maxi:k,
if A[maxi,j]#0 then (
[A[i],A[maxi]]: [A[maxi],A[i]],
A[i]: A[i]/A[i,j],
for u:i+1 thru m do
A[u]: A[u] - A[u,j] * A[i],
i: i+1 ),
j: j+1 ),
A )$
I think this is quite similar to your pseudocode.
Two remarks:
1. Maxima can do parallel assignment, so swapping [a,b]:[b,a] is very easy.
2. If listarith is true (default), Maxima performs basic arithmetrics on lists.
A[i]/A[i,j] divides every list entry by A[i,j].
Now let's see how it works:
(%i14) gaussian_elimination (A) := block( [i, j, m, n, maxi, listarith:true],
...
(%i15) gaussian_elimination (A);
(%o15) matrix([1,-1/3,-1,8/3],[0,1,18/11,-76/11],[0,0,1,-3])
Remarks:
1. The first row differs from (%o7) since echelon doesn't seem to do row swapping.
2. I didn't check if your pseudocode is correct in all cases. I just tried to code it like you
posted it.
3. The above gaussian_elimination changes A by side effects. If you want to avoid this
create a copy of A before processing:
gaussian_elimination (M) := block( [A, i,j, m,n, maxi, listarith:true],
A: copymatrix(M),
...
A )$
HTH
Volker van Nek
Am 10 Aug 2008 um 13:21 hat Angelique Sta maria geschrieben:
>
> How do you program Gaussian Elimination in Maxima?
> I have researched a pseudocode and here it is:
>
> i : 1;
> j : 1;
>
> while (i = m and j = n) do
> Find pivot in column j, starting in row i:
> maxi := i
> for k := i 1 to m do
> if abs(A[k,j]) > abs(A[maxi,j]) then
> maxi := k
> end if
> end for
> if A[maxi,j] ? 0 then
> swap rows i and maxi, but do not change the value of i
> Now A[i,j] will contain the old value of A[maxi,j].
> divide each entry in row i by A[i,j]
> Now A[i,j] will have the value 1.
> for u := i 1 to m do
> subtract A[u,j] * row i from row u
> Now A[u,j] will be 0, since A[u,j] - A[i,j] * A[u,j] = A[u,j] - 1 * A[u,j] = 0.
> end for
> i := i 1
> end if
> j := j 1
> end while
>
> Can you please guide me doing the program?
> Thank you.
>
>
>
> Enrich your blog with Windows Live Writer. Windows Live Writer