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])))$
|
