# Solve a system of time-independent ODE's with vector constants

I have to solve numerically this set of Ordinary Differential Equations

$$\frac{dx_1}{ds} = \frac{1}{x_1} \left[x_2 \left(a + \frac{x_2}{s}\right)-\alpha x_1 z\right]$$

$$\frac{dx_2}{ds} = -\frac{1}{x_1} \left[x_2 \left(b + x_2\right)z\right] - \left(c+\frac{x_2}{s} \right)$$

where $z(s) = \sqrt{x_1^2 +\left(b + x_2 \right)^2}$. The equations are time-independent, $s$ is a distance, $a(s)$, $b(s)$ and $c(s)$ are constant vectors, i.e. $f(s)$; and $\alpha$ is a constant. I've tried to solve using Matlab's ODE45 to no avail because perhaps the fact that the constants are vectors.

Here's the code and some vector examples:

global a b c d
y = ode45(@odeeqns, s, [1 2],[0 0]) % -> I'm not sure about these because I'm working with distance, not time.

function dyds = odeeqns(s,y)
global a b c d
dyds(:,1) = 1./y(1) .*(y(2) .* (a + y(2)./ s) - alpha .* y(1) .* sqrt(y(1).^2 + (b + y(2)).^2));
dyds(:,2) = - (y(2) .* (b + y(2)) .* sqrt(y(1).^2 + (b + y(2)).^2))./y(1) - y(2)./s - c;
plot(dyds);hold on
return


alpha = .0167;
a = [0.0000 0.2985 1.4973 2.4266 2.7838 2.7917 2.6397 2.4295 2.2094 2.0000];
b = [0.0000 0.0298 0.2246 0.4853 0.6960 0.8375 0.9239 0.9718 0.9942 1.0000];
c = [0.0299 0.2615 0.9764 1.4490 1.5680 1.5098 1.3870 1.2499 1.1188 1.0058];
s = [1 2 3 4 5 6 7 8 9 10]; % km


This is the error message:

??? Error using ==> odearguments at 113
YPRIME must return a column vector.

Error in ==> ode45 at 173
[neq, tspan, ntspan, next, t0, tfinal, tdir, y0, f0, odeArgs, odeFcn, ...

Error in ==> windfield_powell at 33
[T,sols] = ode45(@yprime, [1 2],[1 1]);


-
I don't really understand the "type" of your objects. When you say vectors you mean $a(s)$, $b(s)$, $c(s)$ are real numbers? And why do you say that the equations are time independent? ($s$ does appear on the right.) – Sebastien B Feb 11 '13 at 22:31
$s$ is a distance; yes, $a(s).$ are functions (expressed as vectors in Matlab) defined over distance $s$. – Oliver Amundsen Feb 11 '13 at 22:36
Ok, I think I understood your problem. Indeed I think it will be a problem to use ode45. What about coding yourself the explicit euler method or runge-kutta45? Then you can use your vectors. – Sebastien B Feb 12 '13 at 10:55
Do you maybe also have a refenrece to the exercise? Maybe a paper? – k1next Feb 12 '13 at 17:41
nope, did you see the update? – Oliver Amundsen Feb 12 '13 at 17:41

You can try something like that, first define a function

function dx=f(s,x)

a=1.2; b=2.3; c=3.4; alpha=4.5;

dx=zeros(2,1);

z=sqrt(x(1)^2+(b+x(2))^2);

dx(1)=1/x(1)(x(2)(a+x(2)/s)-alpha*x(1)*z);

dx(2)=-1/x(1)(x(2)(b+x(2))*z)-(c+x(2)/s);

and then use ode45 with vectors:

[T,Y] = ode45(@f,[3 12],[1 1])

plot(T,Y(:,1),'-',T,Y(:,2),'-.')

I hope it helps although perhaps you don't want to define the parameters $a$, $b$, $c$, $\alpha$ inside the function $f$

-
The code is helpful, but remember that a, b and c are vectors of constants, not a single constant.... how should I implement that? – Oliver Amundsen Feb 11 '13 at 23:31
@OliverAmundsen can't you implement a function for a? as a = f(s)? – k1next Feb 12 '13 at 11:19
@macydanim, I think that the constant vectors $a$, $b$, and $c$ are already expressed as functions of $s$, in that they're vectors of the same dimensions. That being said, I also tried to integrate the ODE symbolically with Matlab syms, but they're to complicated... any other suggestions? – Oliver Amundsen Feb 12 '13 at 16:02
Just for my comprehension. The ODE system is the one given for $dx_1 /ds$ and $dx_2/ds$ as in the beginng of your question. But actually it isn't $a$ but $a(s)$ and the same with b,c. But these functions are explicitly know? So you could rewrite the ODE system in terms of $x_1,x_2,s$ without using $a(s),b(s),c(s)$ ? Could you give us the functions for $a,b,c$ ? Could you provide your Matlab code? That would help a lot! – k1next Feb 12 '13 at 16:28
Ah maybe I understood something. Do you have prior to running the ODE45 already the values for $a,b,c$ ? Assuming we want to evalute $x_1,x_2$ at $s=0:0.1:1$ then you already have a vector $a$ with $length(a)=10$? And you would like to hand that vector to the function that should be integrate and access the correct values, like for $s=0.2$ you want to access the third element of $a$ as $a(1)$ corresponds to $s=0$ and $a(2)$ to $s=0.1$ ,respectively? Or am I on a complete wrong course? – k1next Feb 12 '13 at 16:35

If you have a function then just add dyds=zeros(size(y))

Looks like this:

function dyds = odeeqns(s,y)

dyds=zeros(size(y))


It should work.

-