Second order system of odes with variable coefficients

[c0der]

Hi,

I have looked everywhere. Can someone please point me in the right direction for solving a system of ODEs with variable coefficients? I managed to solve such system with constant coefficients.

 

[lurflurf]

Solving symbolically or numerically?
Symbolically there is the Peano-Baker Method and use of integrating factors.
it helps if the functions are simple
Have functions in mind?
numerically variable coefficients are no problem if they are well behaved

Have you worked out simple examples like
$$\begin{cases} \mathrm{x}^{\prime}(t)=-t\,\mathrm{y}(t) \\ \mathrm{y}^{\prime}(t)=\phantom{-}t\,\mathrm{x}(t) \end{cases}$$
or
$$\begin{cases} \mathrm{x}^{\prime}(t)=t\,\mathrm{y}(t) \\ \mathrm{y}^{\prime}(t)=\phantom{t}\,\mathrm{x}(t) \end{cases}$$

you can see how quickly things get messy
http://journals.cambridge.org/downl...52a.pdf&code=5491acce57881b30db323feb4fadd1e3

 

[c0der]

Thank you. No I haven't worked out those examples. Do you know of any good examples of ODEs that are solved using the Peano Baker method?

 

[c0der]

Thanks, my equations in the A matrix are:

[ alpha / (b + 2*theta*x) -alpha / (b + 2*theta*x) ]
[ -alpha / (b + 2*theta*(L-x)) alpha / (b + 2*theta*(L-x)) ]

where x is the only variable and this gets complicated enough.

 

[c0der]

After several attempts, I am stuck on the second iteration. I tried integrating in Matlab and an integral does not exist, so I am giving up on this one. Do you think power series would work here e.g. Frobenius method?

 

[lurflurf]

Do you mean
$$
\begin{bmatrix} \\ \mathrm{u}^{\prime}(x) \\ \\ \mathrm{v}^{\prime}(x) \\ \\ \end{bmatrix}
= \alpha\begin{bmatrix}
\phantom{-} \dfrac{1}{b+2\theta x} & -\dfrac{1}{b+2\theta x}
\\
-\dfrac{1}{b+2\theta (L-x)} & \phantom{-}\dfrac{1}{b+2\theta (L-x)}
\end{bmatrix}\begin{bmatrix} \\ \mathrm{u}(x) \\ \\ \mathrm{v}(x) \\ \\ \end{bmatrix}$$

That is an interesting equation
if it is the one you mean it is easily solved
it might help to rename some parameters to make the algebra less messy
write a differential equation for u-v (no u+v) and solve it
write a differential equation for u+v (u-v is now known) and solve it
the final answer has an integral of a radical function that can be left unevaluated or evaluated with hypergeometric function

http://www.wolframalpha.com/input/?.../(b+2theta*(l-t))+alpha*z(t)/(b+2theta*(l-t))

 

[lurflurf]

$$\begin{cases}\frac{\big(\mathrm{u}(x)-\mathrm{v}(x)\big)^\prime}{\mathrm{u}(x)-\mathrm{v}(x)}=\alpha\left(\frac{1}{b+2\, \theta \, x}+\frac{1}{b+2\, \theta \, (L-x)}\right) \\ \frac{\big(\mathrm{u}(x)+\mathrm{v}(x)\big)^\prime}{\mathrm{u}(x)-\mathrm{v}(x)}=\alpha\left(\frac{1}{b+2\, \theta \, x}-\frac{1}{b+2\, \theta \, (L-x)}\right)\end{cases}$$

Hopefully no typos

 

[lurflurf]

Do you know of any good examples of ODEs that are solved using the Peano Baker method?

I do not have a textbook handy.
Besides the above ones some simple examples are
$$\begin{bmatrix} \mathrm{u}^{\prime}(t) \\ \mathrm{v}^{\prime}(t)\end{bmatrix}=
\begin{bmatrix} 0 & 1 \\ a & -a \, t\end{bmatrix}\begin{bmatrix} \mathrm{u}(t) \\ \mathrm{v}(t)\end{bmatrix}$$

$$\begin{bmatrix} \mathrm{u}^{\prime}(t) \\ \mathrm{v}^{\prime}(t)\end{bmatrix}=
\begin{bmatrix} a & e^{-t} \\ 0 & a \end{bmatrix}\begin{bmatrix} \mathrm{u}(t) \\ \mathrm{v}(t)\end{bmatrix}$$

http://arxiv.org/pdf/1011.1775v2.pdf
has two examples

$$\begin{bmatrix} \mathrm{u}^{\prime}(t) \\ \mathrm{v}^{\prime}(t)\end{bmatrix}=
\begin{bmatrix} 0 & t \\ a & 0\end{bmatrix}\begin{bmatrix} \mathrm{u}(t) \\ \mathrm{v}(t)\end{bmatrix}$$

$$\begin{bmatrix} \mathrm{u}^{\prime}(t) \\ \mathrm{v}^{\prime}(t)\end{bmatrix}=
\begin{bmatrix} 1 & t \\ 0 & a \, t\end{bmatrix}\begin{bmatrix} \mathrm{u}(t) \\ \mathrm{v}(t)\end{bmatrix}$$

 

[c0der]

Thanks.

Here is my issue, I even tried solving the equation numerically and the solution blows up incrementally using central differencing.

Here is a simple example system:

d^2y1/dx^2 = alpha * ( y1 - y2 )
d^2y2/dx^2 = alpha * ( y2 - y1 )

As I don't have any initial condition for y(x-h) i.e. y_-1 where h is the timestep size, I am trying to employ the second order forward difference method, but this requires two unknowns y(x+2h) and y(x+h). Any idea what to do?

 

[lurflurf]

What are the initial conditions?
Normally one estimates
y(x+h) and y(x+2h) from the derivative

This one is another where sum difference is helpful

(y1-y2)''=2*alpha * ( y1 - y2 )
(y1+y2)''=0

The solution has e^x so if x is large unless the coefficient is small or zero the solution tends to rapid growth.

 

[c0der]

ICs:

y1(0) = 0
y2(0) = 1000
I tried alpha small 0.0005 and it still grows rapidly, except that this is not the case with the analytical solution

 

[lurflurf]

you need four conditions such as
y1(0), y2(0), y1'(0),and y2'(0)

Probably your numerical method is unstable

 

[c0der]

Of course, I should have included the other 2 for you. y1(0)=0, y2(0)=1000, y1(n)=1000, y2(n)=0. However, this does not seem to work with central differencing no matter what I do as a first derivative BC is required to solve for the first y at -1

 

[lurflurf]

What is n? So this is a boundary value problem so it is a bit different. Some of the simple options include
1)Shooting method
guess y1(1) and y2(1) then you have an initial value problem
compute y1(n) and y2(n) and compare to the given values
repeat until y1(n) and y2(n) match the given values to within tolerance
2)Direct sparse
write 2n-2 equations in 2n-2 variables then direct solve them
2)Iterative sparse
write 2n-2 equations in 2n-2 variables then iterate them

 

[c0der]

To restate the system:

d^2y1/dx^2 = alpha * ( y1 - y2 )
d^2y2/dx^2 = alpha * ( y2 - y1 )

Where alpha is a constant

In numerical form:

[ y1(I+1) - 2y1(I) + y1(I-1) ] / h^2 = alpha * [ (y1(I) - y2(I) ]
[ y2(I+1) - 2y2(I) + y2(I-1) ] / h^2 = alpha * [ (y2(I) - y1(I) ]

Where h is the timestep size

Hence,

y1(I+1) = (2 + alpha * h^2) * y1(I) - y1(I-1) - alpha * h^2 * y2(I)
y2(I+1) = (2 + alpha * h^2) * y2(I) - y2(I-1) - alpha * h^2 * y1(I)

Say I want to use the shooting method, with my BCs y1(0)=0 and y2(0)=sigma0 where sigma0 is any constant. So I need an initial y1'(0) and y2'(0). Given the central difference method [ y1(I+1) - y1(I-1) ] / (2h) where h is the timestep. Thus initially, I have:

a) [ y1(1) - y1(-1) ] / 2h = a1 for my initial guess for the derivative, hence y1(-1) = y1(1) - 2*h*a1
b) Similarly, y2(-1) = y2(1) - 2*h*a2
c) When I=0, substituting from (a) and (b) into the above numerical expressions:

y1(1) = (2 + alpha * h^2) * y1(0) - y1(-1) - alpha * h^2 * y2(0)
y1(1) + y1(1) - 2*h*a1 = (2 + alpha * h^2) * y1(0) - alpha * h^2 * y2(0)
y1(1) = [ (2 + alpha * h^2) * y1(0) - alpha * h^2 * y2(0) - 2 * h * a1 ] / 2

y2(1) = (2 + alpha * h^2) * y2(0) - y2(-1) - alpha * h^2 * y1(0)
y2(1) + y2(1) - 2*h*a2 = (2 + alpha * h^2) * y2(0) - alpha * h^2 * y1(0)
y2(1) = [ (2 + alpha * h^2) * y2(0) - alpha * h^2 * y1(0) - 2 * h * a2 ] / 2

For each timestep at I=1 to N

% Iterate here until difference between IVP and BVP is small

y1(I+1) = (2 + alpha * h^2) * y1(I) - y1(I-1) - alpha * h^2 * y2(I)
y2(I+1) = (2 + alpha * h^2) * y2(I) - y2(I-1) - alpha * h^2 * y1(I)

End for

 

[the_wolfman]

One of the terms in the analytic solution is exponentially growing. This term could be giving you problems. In order to determine if this is the problem we need to know both the length of the system and your mesh spacing. You want the mesh spacing to be smaller than the characteristic distance of the exponential term. You also don't want the domain size to be much larger than this characteristic distance.

Also what values are you using for as your initial guesses for $y_1'$ and $y_2'$?

It may be worth wild to further simplify the system.
What happens if you set $\alpha=0$. Is you system stable? Do reproduce the analytic solution for this case?

Second, what happens if you try solving the equation
$y_1''=\alpha y_1$?

 

[c0der]

Thanks for the reply. Does reduction of order work on systems of ODEs with variable coefficients?

 

[lurflurf]

^Yes reduction of order works on systems of ODEs with variable coefficients.

 

[c0der]

I thought so. Here is what I did and the solution didn't pan out.


d^2(sigma1)/dx^2 = -(2*theta) / (b + 2*theta*x) * d(sigma1)/dx + alpha / (b + 2*theta*x) * (sigma1 - sigma2)
d^2(sigma2)/dx^2 = (2*theta) / [ b + 2*theta*(L-x) ] * d(sigma2)/dx + alpha / [ b + 2*theta*(L-x) ] * (sigma2 - sigma1)


1) Reducing these to a first order system:



Let y1 = sigma1, y2 = d(sigma1)/dx

d(y1)/dx = d(digma1)/dx = y2 and;

d(y2)/dx = d^2(sigma1)/dx^2 = -(2*theta) / (b + 2*theta*x) * d(sigma1)/dx + alpha / (b + 2*theta*x) * (sigma1 - sigma2)



Let y3 = sigma2, y4 = d(sigma2)/dx

d(y3)/dx = d(sigma2)/dx = y4 and;

d(y4)/dx = d^2(sigma2)/dx^2 = (2*theta) / [ b + 2*theta*(L-x) ] * d(sigma2)/dx + alpha / [ b + 2*theta*(L-x) ] * (sigma2 - sigma1)



Then the first order system is:



d(y1)/dx = y2

d(y2)/dx = -(2*theta) / (b + 2*theta*x) * y2 + alpha / (b + 2*theta*x) * (y1 - y3)

d(y3)/dx = y4

d(y4)/dx = (2*theta) / [ b + 2*theta*(L-x) ] * y4 + alpha / [ b + 2*theta*(L-x) ] * (y3 - y1)



I then take the Laplace transforms, applying the BCs (y1(0) = sigma1(0) = 0, y3(0) = sigma2(0) = sigma0) to get:



s*Y1(s) - Y2(s) = 0

-alpha*Y1(s) + b*s*Y2(s) + alpha*Y3(s) = b*y2(0)

s*Y3(s) - Y4(s) = sigma0

alpha*Y1(s) - alpha*Y3(s) + (bs + 2*theta*L*s)*Y4(s) = (b + 2*theta*L)*y4(0)



Then I solve the system using Gauss elimination and back substitute to make sure I got the correct answer, which I did.



Finally, I take the inverse Laplace transform in Matlab and got a solution for y1 and y3 which represent sigma1 and sigma2 respectively. Differentiating them appropriately and plugging them back into the original ODEs does not get me the correct answer.



Even though my boundary conditions for the first derivative are unknown at x=0, I would still expect that these steps yield a correct solution because the unknowns at x=0 are still constants y2(0) and y4(0).



I would appreciate any suggestions.