Looking for an elegant way to solve a (system of) ODEs/functional equations with undetermined coefficients

Question

I want to solve an ODE using undetermined coefficients/guess-and-verify and am looking for an elegant way to use this technique. I am running into a few ugly warts. Here is a basic implementation for a simple example:

ode = r v[z] == Exp[z] + g v'[z] + s^2/2 v''[z]
guess = {v[z_] -> a Exp[z] + b Exp[c z]} (* for some a,b,c*)

(* Inelegance 1: how can I substitute the guess into the ode directly?
For example, ode //. guess doesn't work  *)
guessderivs = 
  Union[guess, {v'[z_] -> (D[v[z] /. guess, z]), v''[z_] -> (D[v[z] /. guess, z, z])}]
(* ugly, but complete now.  Substitute*)

odesub = ode //. guessderivs // FullSimplify
(* => e^z (1 + a (g - r)) + b e^(c z) (c g - r) == 0 *)

(*Inelegance #2: This is an easy example, but is there a standard way to
get equations to solve for all z?  SolveAlways[odesub,z] didn't work *)

(* Doing it manually works, of course*)
    a$sub = First@Solve[(2 + a (2 g - 2 r + s^2)) == 0, a]
    c$sub = Solve[(2 c g - 2 r + c^2 s^2) == 0,c]

(* The last constant, b, would come from an initial value. mine is v'[0]==0 for example *)

To summarize where this is ugly:

1. How to make a guess substitution of a functional form, and have the derivatives work (i.e. v'[z] //. {v[z_]-> a Exp[b z]} but actually substitute the derivatives) I recognize the problem here that FullForm[v'[z]] == Derivative[1][v][z] so the rule doesn't match, but is there a nice way to replace these generally?

2. How to get the system of equations (in general) with undetermined coefficients? I thought that SolveAlways was supposed to do this?

EDIT: It looks like the cleanest method is from @RiemannZeta which is to give a list of linearly independent basis vectors (anonymous as described by @belisarius), substitute, and find the coefficients.

For this setup, let me give the following goal for some generic UndeterminedCoefficientsSystem function:

ode = {r v[z] == Exp[z] + g v'[z] + s^2/2 v''[z]};
ic = {v'[0] == 0};
equations = Union[ode, ic]; (* In reality, I will be trying a system of equations *)
v$basis = {Exp[#]&, Exp[d #]&, # &}; (* with an addition to test, which should be 0 here *)
    coefficientsystem = UndeterminedCoefficientSystem[{{v, v$basis}}, equations};

(* Assume that this function generates a list of coefficients of the form C[1][1], C[1][2],... for the entire basis then it substitutes into the set of equations  *)
coefficientsystem == {C[1][1] == -1 - g C[1][1] + r C[1][1] - 1/2 s^2 C[1][1],
C[1][2] == -d g C[1][2] + r C[1][2] - 1/2 d^2 s^2 C[1][2],
C[1][3] == r C[1][3],
C[1][1] + d C[1][2] + C[1][3] == 0
} (* The last equation comes from the v'[0]==0 *)

So basically, the function does the following:

1. For arguments, takes a list of function name with a list of basis as anonymous functions. Here it is only 1, but I need to use this to solve a more complicated system of ODEs. Index the functions by j.
2. Take each function with index j and basis of length i of anonymous functions. Generate a list of coefficients for C[i][j] of these.
3. For each function and basis, create a guess anonymous function as the linear combination of these, with the C[i][j] coefficients. Create a substitution rule given the function name to the linear combination of basis.
4. For each given equation, substitute in the concatenated list of guesses (for each j)
5. For each basis, for each equation, find the coefficient with @RiemannZeta technique and append the whole list (if there are different basis for the different functions, it is reasonable to just look for the coefficients of the union of all anoymous basis)
6. Return the whole list of equations.

Does that make sense as a generic function? It should work for any sort of functional equations, and not just ODEs.

(I am getting stuck very early here. For example, we need to have anonymous guesses (to ensure that we can both take derivatives and evaluate at initial values such as v'[0] == 0) but how do I take two anonymous functions and concatenate a linear combination to form a new one?)

Answer 1

To see the steps you can just ask Wolfram|Alpha :)

WolframAlpha["r v[z] == Exp[z] + g v'[z] + s^2/2 v''[z]"]

But if you're wanting to do it yourself, here's what I recommend.

1. I'd do the following:

   guess = a Exp[z] + b Exp[c z];
   odesub = ode /. {v[z] -> guess, Derivative[n_][v][z] :> D[guess, {z,n}]} // FullSimplify
   
   (* b E^(c z) (2 c g - 2 r + c^2 s^2) + E^z (2 + a (2 g - 2 r + s^2)) == 0 *)
   

2. For this step, I think you want to equate coefficients of your linearly independent basis functions to zero.

   coeffs = Coefficient[Subtract @@ odesub, {Exp[z], Exp[c z]}]
   (* {2 + 2 a g - 2 a r + a s^2, 2 b c g - 2 b r + b c^2 s^2} *)
   
   Solve[Thread[coeffs == 0], {a, c}]
   (* {{a -> -2/(2 g - 2 r + s^2), c -> (-g - Sqrt[g^2 + 2 r s^2])/s^2},
       {a -> -2/(2 g - 2 r + s^2), c -> (-g + Sqrt[g^2 + 2 r s^2])/s^2}} *)
   

Answer 2

ode = r v[z] == Exp[z] + g v'[z] + s^2/2 v''[z];
guess = {v -> (a Exp[#] + b Exp[c #] &)};
Resolve[ForAll[z, ode /. guess], Reals]
(*
  ...a huge mess, after a while ...
*)