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

Question

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

---

Any help please? Thanks!

Answer 1

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$