I came up with my own solution to analytically solving very simple 1st and
second order ODEs (i.e. returns the function/ verb solution) with boundary
conditions.
For example to solve y'' +2y' +y = 0, with y(0) = 1 and y'(0) = 0
I do
(0 1; 0 0; 1 0) SOLVE_ODE 1 2 1
0.999999999999999667&*@:(^@:(_1x&*)) + 0.999999999999999667&*@:(] * ^@:(_1x&*))
NB. i wish it would round to one
It also took care of the repeated root in the second summand.
The three elements of the boxed array in the left argument of SOLVE_ODE are:
the derivatives of the boundary conditions; the arguments of the boundary
conditions; the results of the functions on the arguments.
Another example:
y''+4y'-3y=0, with y(1) = 5 and y''(-1) = 5
solution =: (0 2; 1 _1; 5 5) SOLVE_ODE _3 4 1
Result is:
0.00196954874057786878&*@:(^@:(_4.64575131106459072&*)) +
2.62133261178797961&*@:(^@:(0.645751311064590605&*))
solution 1
gives 5.
My solution also handles some first order linear ODEs :
e.g. 2y' +y = 0, with y'(0) = 1
a=: ( 1; 0 ; 1 ) SOLVE_ODE 1 2
returns
_2&(^@:(_0.5&*))
My solution didn't follow the (better) matrix approach. I decided to solve it
this way just to see how to do such a thing in J. It is also clumsy handling
complex solutions, which in theory should be able to be reduced to sin / cosine
functions. But I haven't figured a way to do that.
It also cannot handle non-zero arguments on the right side of the ODE.
e.g. y'' + y' - y = cos(x)
is out of the question.
Anyway, I have added my source below. I would appreciate any comments for
improving the general J-ness of my code. If you look at the end of
BOUND_CONDITIONS, I ran into a problem trying to return the verb on the last
line (before the final end.) Returning conditions&r1 did not give me the actual
verb solution. Any explanation of what my mistake was would be appreciated.
Also, I was able to write a shorter program by calculating both summands
together (i.e. both exponentials). The problem then is that I cannot call d. on
this
e.g.
F =: +/@:((1 2)&*)@:^@:(_1 1&*)
(F d. 1) 1 NB. differentiate and evaluate at x = 1.
This gives an error. It seems d. only handles scalar arguments of verbs (?).
And I really want to be able to differentiate my solution.
My code:
ROOTS=: >@:(1&{)@p.
G=: 1 : '^@:(m&*)'
SOLVE=: 1 : ' (ROOTS m) G2'
SOLVEn=: 2 : '(n{ (ROOTS m)) G2'
NB. Gets the solution
G2=: 1 : 0
sols=: m
if. (# sols ) = 2 do.
if. =/ sols do.
((>:)*(^@:(_1x&*)))
elseif. 1 do. ^@:(m&*)
end.
else.
^@:(m&*)
end.
)
NB. Returns general solution verb, without constant coefficients.
SOLVE=: 1 : ' (ROOTS m) G2'
NB. Solution for duplicate roots.
SOLVE_DUPLICATE=: 2 : 0
root=. 0{ (ROOTS n)
((0{m)&*@:(^@:(root&*)))+ ((1{m)&*@:(]*(^@:(root&*))))
)
NB. Boundary Conditions. Returns the coefficients
NB. of the soluton's summands, calculated from
NB. the given boundary conditions.
BOUND_CONDITIONS=: conjunction define
coeffs=. y NB. polynomial coefficients
deriv=. >0{ x NB. the orders of derivatives
val=. >1{ x NB. values to put in
res=. >2{ x NB. the values' results, i.e. y^(deriv)(val) = res
if. (# coeffs) = 3 do. NB. Quadratic ODE (e.g. y''+y'+y = 0)
if. =/ (ROOTS y) do.
c1=. 1 :'^@:((0{ (ROOTS m))&*)'
c2=. 1 :']*(^@:((0{ (ROOTS m))&*))'
r1=. y c1
r2=. y c2
NB. non-repeating roots. Get the two summands of the
NB. solution independently.
else.
r1=. y SOLVEn 0
r2=. y SOLVEn 1
end.
NB. solutions - for quadratics, differentiate the
NB. summands and input the values.
NB. This gets four numbers, we can put into a matrix to find
NB. the 2 coefficients of the summands.
dr1=. (r1 d. (0{deriv)) (0{val)
dr2=. (r2 d. (0{deriv)) (0{val)
dr3=. (r1 d. (1{deriv)) (1{val)
dr4=. (r2 d. (1{deriv)) (1{val)
NB. matrixify, get the coefficient constants
mat=. (2 2) $ dr1, dr2, dr3, dr4
mat=. %. mat
constants=. mat (+/ . *) res
constants NB. these are the two constants for summands 1 and 2.
if. =/ (ROOTS y) do.
((0{constants)&*@:(y c1)) + ((1{constants)&*@:(y c2))
else.
((0{constants)&*@:(y SOLVEn 0)) + ((1{constants)&*@:(y SOLVEn 1))
end.
elseif. (# coeffs) = 2 do. NB. order 1 eqn (e.g. y'+y = 0)
r1=. y SOLVEn 0 NB. ge the general solution, without coefficient.
dr1=. (r1 d. (0{deriv)) (0{val)
constant=. res % dr1 NB. coefficient.
constant &(y SOLVEn 0) NB. returning constant&r1 doesn't work, so (y SOLVEn
0) needed to be called again here.
end.
)
SOLVE_ODE =: conjunction define
if. (# n) = 3 do.
if. =/( ROOTS n) do.
(m BOUND_CONDITIONS n)
else.
(m BOUND_CONDITIONS n )
end.
elseif. (# n ) = 2 do.
(m BOUND_CONDITIONS n )
end.
)
> Date: Sun, 15 Mar 2015 21:14:17 -0500
> From: [email protected]
> To: [email protected]
> Subject: Re: [Jprogramming] Solving Differential Eqns with J
>
> Use an adverb:
>
> G =: 1 : '^@:(m&*)'
> _2 1 G
> ^@:(_2 1&*)
>
> --Kip Murray
>
> On Sunday, March 15, 2015, Jon Hough <[email protected]> wrote:
>
> > Kip, Raul,Thanks for replying.
> > The matrix way would probably be better.
> > As a follow up question, I am having trouble creating an
> > adverb/conjunction to build solutions from.
> > if I define the conjunction
> >
> > F =: 2 : '^@:(n&*)'
> >
> >
> > and I give, as my noun argument, the solutions to some polynomial (e.g.
> > >1{ p. _2 1 1)
> >
> > Then I need to do
> >
> > '' F _2 1
> >
> >
> >
> >
> > ^@:(_2 1&*) NB. the two solutions (without constants)
> >
> >
> >
> >
> > Which gives me my verb "solution". Please note that I had to give '' as
> > the left argument to F to return my verb.
> >
> >
> > Is there some other way of doing this?
> > i.e. I want some "part of speech" (conjunction or adverb I suppose) that
> > takes
> > a noun and returns a verb.
> >
> >
> > It seems conjunctions take verbs and nouns and return verbs
> >
> >
> > Conj: V x N --> V
> >
> >
> > but since, in my case for F, the verb is superfluous, is there any "part
> > of speech" that takes a noun and returns a verb?
> >
> >
> > I hope my question made sense. My motivation is, I don't see why I should
> > need to pass an empty verb '' to F to
> > return a verb. I only want to give F my noun and from that, I want to get
> > a verb.
> >
> >
> > Regards,
> > Jon
> >
> >
> >
> > > From: [email protected] <javascript:;>
> > > Date: Sun, 15 Mar 2015 14:46:35 -0400
> > > To: [email protected] <javascript:;>
> > > Subject: Re: [Jprogramming] Solving Differential Eqns with J
> > >
> > > Oops, I should have seen that.
> > >
> > > Thanks,
> > >
> > > --
> > > Raul
> > >
> > > On Sun, Mar 15, 2015 at 2:45 PM, Kip Murray <[email protected]
> > <javascript:;>> wrote:
> > > > Raul, where you have _1 * y you need _2 * y
> > > >
> > > > --Kip
> > > >
> > > > On Sunday, March 15, 2015, Raul Miller <[email protected]
> > <javascript:;>> wrote:
> > > >
> > > >> On Sat, Mar 14, 2015 at 7:11 AM, Kip Murray <[email protected]
> > <javascript:;>
> > > >> <javascript:;>> wrote:
> > > >> > I recommend a matrix approach involving eigenvalues and
> > eigenvectors,
> > > >> > and assume you have a way to find eigenvalues and eigenvectors for
> > an n
> > > >> by
> > > >> > n matrix which has n linearly independent eigenvectors. I confine
> > my
> > > >> > attention to linear differential equations with constant
> > coefficients and
> > > >> > right hand side 0, for example
> > > >> >
> > > >> > u''' - u'' - 4 u' + 4 u = 0 with initial conditions u(0) = 2 ,
> > u'(0) =
> > > >> _1
> > > >> > , and u''(0) = 5
> > > >> >
> > > >> > This has the solution u(t) = ( ^ _2 * t ) + ( ^ t ) where I use
> > a mix
> > > >> of
> > > >> > conventional and J notation.
> > > >>
> > > >> I decided to try converting your notes here to J, and I ran into a
> > snag.
> > > >>
> > > >> First, I tried expressing your initial constraints in J:
> > > >>
> > > >> constraint=: 1 :0
> > > >> if. 0 = y do.
> > > >> assert. 2 = u y
> > > >> assert. _1 = u D. 1 y
> > > >> assert. 5 = u D. 2 y
> > > >> end.
> > > >> assert. 0= (u D.3 + (_1*u D.2) + (_4 * u D. 1) + 4 * u) y
> > > >> )
> > > >>
> > > >> And then I expressed your solution in J:
> > > >>
> > > >> solution=:3 :0
> > > >> (^ _1 * y) + ^ y
> > > >> )
> > > >>
> > > >> Or, alternatively:
> > > >>
> > > >> tacitsolution=: +&^ -
> > > >>
> > > >> And then I tested my code to see if I had gotten it right:
> > > >>
> > > >> solution constraint 0
> > > >> |assertion failure
> > > >> | _1=u D.1 y
> > > >>
> > > >> Looking at this:
> > > >> solution D.1 (0)
> > > >> 9.76996e_8
> > > >> tacitsolution D.1 (0)
> > > >> 9.76996e_8
> > > >>
> > > >> working this through manually,
> > > >>
> > > >> dsolution=: 3 :0
> > > >> (_1 * ^ _1*y) + ^ y
> > > >> )
> > > >>
> > > >> dsolution 0
> > > >> 0
> > > >>
> > > >> So J's implementation of D. isn't perfect, but I'm wondering if you
> > > >> might not also have a typo somewhere in your presentation here?
> > > >>
> > > >> Thanks,
> > > >>
> > > >> --
> > > >> Raul
> > > >> ----------------------------------------------------------------------
> > > >> For information about J forums see
> > http://www.jsoftware.com/forums.htm
> > > >>
> > > >
> > > >
> > > > --
> > > > Sent from Gmail Mobile
> > > > ----------------------------------------------------------------------
> > > > For information about J forums see http://www.jsoftware.com/forums.htm
> > > ----------------------------------------------------------------------
> > > For information about J forums see http://www.jsoftware.com/forums.htm
> >
> > ----------------------------------------------------------------------
> > For information about J forums see http://www.jsoftware.com/forums.htm
> >
>
>
> --
> Sent from Gmail Mobile
> ----------------------------------------------------------------------
> For information about J forums see http://www.jsoftware.com/forums.htm
----------------------------------------------------------------------
For information about J forums see http://www.jsoftware.com/forums.htm
