1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
|
/* THIS LITTLE PACKAGE SOLVES FIRST ORDER LINEAR
ORDINARY DIFFERENTIAL EQUATIONS SUBJECT TO A
BOUNDARY CONDITION (B.C.) AT AN INITIAL POINT
THE CALLING PROCEDURE IS
IVPSOL(DIFFEQ,Y,X,A,BCEQ);
WHERE DIFFEQ IS THE DIFFERENTIAL EQUATION TO BE SOLVED
(E.G.: X*'DIFF(Y,X)+2*Y=1)
Y IS THE DEPENDENT VARIABLE
X IS THE INDEPENDENT VARIABLE
A IS THE POINT AT WHICH THE B.C. IS APPLIED
BCEQ IS THE B.C. EQUATION (E.G.: Y=2)
*/
SOLDIFF(EQ1,Y,X,A):=
BLOCK([EQ2,A1,A2,A3,A4,A5],
EQ2: LHS(EQ1)-RHS(EQ1),
A1: RATCOEF(EQ2,'DIFF(Y,X)),
A2: RATCOEF(EQ2,Y),
A3: EQ2-A1*'DIFF(Y,X)-A2*Y,
A4: INTEGRATE(A2/A1,X),
A5: INTEGRATE(%E**A4*A3/A1,X),
'CONST1/%E**A4-A5/%E**A4)$
IVPSOL(DIFFEQ,Y,X,A,BCEQ):=
BLOCK([Z,DERY,DYA,YA,CO1],
Z: SOLDIFF(DIFFEQ,Y,X,A),
DERY:DIFF(Z,X,1),
DYA: SUBST(A,X,DERY),
YA: SUBST(A,X,Z),
CO1: SOLVE(SUBST(YA,Y,SUBST(DYA,'DIFF(Y,X),BCEQ)),'CONST1),
IF length(co1)>1 /* LISTP(CO1) */ THEN
[PRINT("SOLUTIONS ARE NOT UNIQUE, POSSIBLE SOLUTIONS ARE:
"),
MAP(LAMBDA([L1],apply('DISPLAY,[L1])),CO1), RETURN(Z)]
ELSE if co1=[] then return(print("No solution for given initial
condition"),z) else SUBST(RHS(CO1[1]),'CONST1,Z))$
/* FIND RIGHT EIGENVECTORS WITH OTHER COMPONENT = 1
I.E. X IS SOLUTION OF M*X=MU*X */
EVEC(M,MU,MODES):=block([equations,unknowns,x],
(FOR L:1 THRU MODES DO[
EQUATIONS:[],
UNKNOWNS:[],
/* KILL(X), */
X[L]:1,
FOR I:1 THRU MODES DO IF NOT (I=L) THEN [
UNKNOWNS:CONS(X[I],UNKNOWNS),
EQUATIONS:CONS(SUM(M[I,J]*X[J],J,1,MODES)
=MU[L]*X[I],EQUATIONS)],
PRINT(SOLVE(EQUATIONS,UNKNOWNS))]))$
/* CONSTRUCT 3 X 3 MATRIX WITH DESIRED EIGENVALUES */
CONSMAT3(L1,L2,L3):=BLOCK([A,B,MAT3,CP,CPM,SOL1,L],
MAT3:MATRIX([1,2,3],[B,3,A],[1,1,L1+L2+L3-4]),
CP:EXPAND(-(L-L1)*(L-L2)*(L-L3)),
CPM:CHARPOLY(MAT3,L),GLOBALSOLVE:TRUE,
SOL1:SOLVE([RATCOEF(CP,L)=RATCOEF(CPM,L),
RATCOEF(CP,L,0)=RATCOEF(CPM,L,0)],[A,B]),
RETURN(apply('EV,[MAT3])))$
|
