Discover MakerZone

MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

Learn more

Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

Apply Today

Thread Subject:
how to Solve coupled ODEs

Subject: how to Solve coupled ODEs

From: Henry

Date: 6 Feb, 2013 21:59:09

Message: 1 of 14

Hi there,

I'm having a hard time solving a set of 2 coupled ODEs:

dx/ds = 1/x * [y* (d + y/s) - a*x*f]
dy/ds = 1/x * [-y * (b + y) * f] - y/s - c

where f = sqrt(x^2 + y^2*(e + y)), and all constants are vectors of (100 x 1). The independent variable is "s".

%----------------------------------------------------------------------
global a b c d e
y = ode45(@odeeqns, s, [1 1])

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

Any help please?
Thanks!

Subject: how to Solve coupled ODEs

From: Torsten

Date: 7 Feb, 2013 08:15:23

Message: 2 of 14

"Henry" wrote in message <keujnd$h64$1@newscl01ah.mathworks.com>...
> Hi there,
>
> I'm having a hard time solving a set of 2 coupled ODEs:
>
> dx/ds = 1/x * [y* (d + y/s) - a*x*f]
> dy/ds = 1/x * [-y * (b + y) * f] - y/s - c
>
> where f = sqrt(x^2 + y^2*(e + y)), and all constants are vectors of (100 x 1). The independent variable is "s".
>
> %----------------------------------------------------------------------
> global a b c d e
> y = ode45(@odeeqns, s, [1 1])
>
> function dyds = odeeqns(s,y)
> global a b c d e
> dyds(:,1) = 1./y(1) .*(y(2) .* (d + s ./ y(2)) - a .* y(1) .* sqrt(y(1).^2 + (e + y(2)).^2));
> dyds(:,2) = - (y(2) .* (b + y(2)) .* sqrt(y(1).^2 + (e + y(2)).^2))./y(1) - y(2)./s - c;
> return
> %----------------------------------------------------------------------
>
> Any help please?
> Thanks!

Take a look at example 1 under
http://www.mathworks.de/de/help/matlab/ref/ode23.html
and you'll recognize the two (or three) errors you made in your small piece of code
from above.

Best wishes
Torsten.

Subject: how to Solve coupled ODEs

From: Henry

Date: 7 Feb, 2013 17:11:08

Message: 3 of 14

Hi Torsten,
Would you kindly indicate which are the mistakes in my code?
thanks! Henry

Subject: how to Solve coupled ODEs

From: Torsten

Date: 8 Feb, 2013 07:47:07

Message: 4 of 14

"Henry" wrote in message <kf0n7c$k5i$1@newscl01ah.mathworks.com>...
> Hi Torsten,
> Would you kindly indicate which are the mistakes in my code?
> thanks! Henry

Just compare line for line:

 global a b c d e
 a=...;
 b=...;
 c=...;
 d=...;
 e=...;
 [t,y] = ode45(@odeeqns, [0 3], [1 1]);
 
 function dyds = odeeqns(s,y)
 global a b c d e
 dyds = zeros(2,1)
 dyds(1,1) = 1/y(1) *(y(2)*(d+s /y(2))-a *y(1)*sqrt(y(1)^2+(e+y(2))^2));
 dyds(2,1) = -(y(2)*(b+y(2))*sqrt(y(1)^2+(e+y(2))^2))/y(1)-y(2)/s - c;

Best wishes
Torsten.

Subject: (yet unsolved) how to Solve coupled ODEs

From: Henry

Date: 11 Feb, 2013 17:10:27

Message: 5 of 14

Problem not yet solves, any helps please?

Subject: (yet unsolved) how to Solve coupled ODEs

From: Torsten

Date: 12 Feb, 2013 16:07:18

Message: 6 of 14

"Henry" wrote in message <kfb8m3$mom$1@newscl01ah.mathworks.com>...
> Problem not yet solves, any helps please?

Just compare line for line:

 global a b c d e
 a=...;
 b=...;
 c=...;
 d=...;
 e=...;
 [t,y] = ode45(@odeeqns, [0 3], [1 1]);
 
 function dyds = odeeqns(s,y)
 global a b c d e
 dyds = zeros(2,1)
 dyds(1,1) = 1/y(1) *(y(2)*(d+s /y(2))-a *y(1)*sqrt(y(1)^2+(e+y(2))^2));
 dyds(2,1) = -(y(2)*(b+y(2))*sqrt(y(1)^2+(e+y(2))^2))/y(1)-y(2)/s - c;

Best wishes
Torsten.

If you start with s=0 as I suggested in the call to ode45, you will get problems in the evaluation of dyds(2,1) because of the term y(2)/s (division by zero).
Please supply the values of the parameters a,b,c,d and e, the interval of integration and the error message you received.

Best wishes
Torsten.

Subject: (yet unsolved) how to Solve coupled ODEs

From: Henry

Date: 12 Feb, 2013 17:17:31

Message: 7 of 14

Hi Torsten, here it is. I updated the function dyds removing a couple of typos, and provide the vectors for only 10 elements (s in KM).

-----
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) .* (d + y(2)./ s) - a .* 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
----
a = .0167;
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];
d = [0.0000 0.2985 1.4973 2.4266 2.7838 2.7917 2.6397 2.4295 2.2094 2.0000];
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]);
-----
Thanks!

Subject: (yet unsolved) how to Solve coupled ODEs

From: Steven_Lord

Date: 12 Feb, 2013 18:07:07

Message: 8 of 14



"Henry " <mflacius1521@gmail.com> wrote in message
news:kfdtfb$719$1@newscl01ah.mathworks.com...
> Hi Torsten, here it is. I updated the function dyds removing a couple of
> typos, and provide the vectors for only 10 elements (s in KM).
>
> -----
> 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) .* (d + y(2)./ s) - a .* 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
> ----
> a = .0167;
> 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];
> d = [0.0000 0.2985 1.4973 2.4266 2.7838 2.7917 2.6397 2.4295 2.2094
> 2.0000];
> 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.

Well, DOES your function return a column vector? As written above it does
not. Modify it so it does and try again. [You will likely run into a
different error, because you haven't provided as many initial conditions as
you have ODEs. You'll need to determine the correct initial conditions on
your own since you haven't provided the group with a mathematical or
physical statement of the problem underlying this system of ODEs.]

--
Steve Lord
slord@mathworks.com
To contact Technical Support use the Contact Us link on
http://www.mathworks.com

Subject: (yet unsolved) how to Solve coupled ODEs

From: Henry

Date: 12 Feb, 2013 18:26:05

Message: 9 of 14

hi, I tried the following

[T,sols] = ode45(@yprime, [1 20],0.000001*ones(1,40)');

function dyds = odeeqns(s,y)
   global a b c d
   dyds1 = 1./y(1) .*(y(2) .* (d + y(2)./ s) - a .* y(1) .* sqrt(y(1).^2 + (b + y(2)).^2));
   dyds2 = - (y(2) .* (b + y(2)) .* sqrt(y(1).^2 + (b + y(2)).^2))./y(1) - y(2)./s - c;
   dyds = [dyds1 dyds2];
size(dyds) % 1 x 40
   plot(dyds);hold on
return

same error message... The funny thing is that the code plots the vectors... but then shows the error message...

the problem is a profile described over a distance from 1 to 20 kms in this case...What am I doing wrong?

Subject: (yet unsolved) how to Solve coupled ODEs

From: Steven_Lord

Date: 12 Feb, 2013 20:04:40

Message: 10 of 14



"Henry " <mflacius1521@gmail.com> wrote in message
news:kfe1ft$mlu$1@newscl01ah.mathworks.com...
> hi, I tried the following
>
> [T,sols] = ode45(@yprime, [1 20],0.000001*ones(1,40)');
>
> function dyds = odeeqns(s,y)
> global a b c d
> dyds1 = 1./y(1) .*(y(2) .* (d + y(2)./ s) - a .* y(1) .* sqrt(y(1).^2 +
> (b + y(2)).^2));
> dyds2 = - (y(2) .* (b + y(2)) .* sqrt(y(1).^2 + (b + y(2)).^2))./y(1) -
> y(2)./s - c;
> dyds = [dyds1 dyds2];
> size(dyds) % 1 x 40

An array with 1 row and 40 columns is a vector, but it is not a COLUMN
vector.

An array with 40 rows and 1 column is a column vector.

--
Steve Lord
slord@mathworks.com
To contact Technical Support use the Contact Us link on
http://www.mathworks.com

Subject: (yet unsolved) how to Solve coupled ODEs

From: Henry

Date: 12 Feb, 2013 21:08:08

Message: 11 of 14

I modified the code to:

tspan = linspace(1, 2,100);
y = 0.0001*ones(1,200);
[T,sols] = ode45(@yprime, tspan,y);
and
dyds = [dyds1' ; dyds2'];

and it now runs, thanks.

But I have this question: My odes are time-independent, function of "s" which is a distance, but since the ode45 will not work without a tspan, I came up with an arbitrary tspan. It looks as if ode45 interprets that "s" is a time? I'm confused as to how interpret the result of the routine which has dimensions of 100x200. I was expecting 1 x 200, i.e. independent of time.
thanks

Subject: (yet unsolved) how to Solve coupled ODEs

From: Torsten

Date: 13 Feb, 2013 08:01:09

Message: 12 of 14

"Henry" wrote in message <kfeavo$2f$1@newscl01ah.mathworks.com>...
> I modified the code to:
>
> tspan = linspace(1, 2,100);
> y = 0.0001*ones(1,200);
> [T,sols] = ode45(@yprime, tspan,y);
> and
> dyds = [dyds1' ; dyds2'];
>
> and it now runs, thanks.
>
> But I have this question: My odes are time-independent, function of "s" which is a distance, but since the ode45 will not work without a tspan, I came up with an arbitrary tspan. It looks as if ode45 interprets that "s" is a time? I'm confused as to how interpret the result of the routine which has dimensions of 100x200. I was expecting 1 x 200, i.e. independent of time.
> thanks

Just interpret your results for y(1) and y(2) as functions of distance, not of time.
So "tspan" in your case is a distance vector, not a time vector.

Best wishes
Torsten.

Subject: (yet unsolved) how to Solve coupled ODEs

From: Torsten

Date: 13 Feb, 2013 09:11:05

Message: 13 of 14

"Henry" wrote in message <kfeavo$2f$1@newscl01ah.mathworks.com>...
> I modified the code to:
>
> tspan = linspace(1, 2,100);
> y = 0.0001*ones(1,200);
> [T,sols] = ode45(@yprime, tspan,y);
> and
> dyds = [dyds1' ; dyds2'];
>
> and it now runs, thanks.
>
> But I have this question: My odes are time-independent, function of "s" which is a distance, but since the ode45 will not work without a tspan, I came up with an arbitrary tspan. It looks as if ode45 interprets that "s" is a time? I'm confused as to how interpret the result of the routine which has dimensions of 100x200. I was expecting 1 x 200, i.e. independent of time.
> thanks

By the way:
I suspect that you want to set different values of a,b,c and d for different values of s and solve for only _two_ solutions y(1) and y(2). That's not what you do in your code above.
You will have to use MATLAB's interp1 function to get values for a,b,c and d at the distance s where the solver needs the values of the derivatives dyds(1) and dyds(2).

Or do you want to make a parameter study, i.e. different solutions for different
combinations of the parameters a,b,c and d ? Then your code above
is also incorrect because you always refer to y(1) and y(2) in the evaluation of the derivatives, thus the solution for only the _first_ parameter set.

We need more information about your underlying problem in order to give
useful advice.

Best wishes
Torsten.

Subject: (yet unsolved) how to Solve coupled ODEs

From: Henry

Date: 13 Feb, 2013 20:41:14

Message: 14 of 14

Hi Torsten,

As you mention in your second suggestion, I guess I need to get "different solutions for different combinations of the parameters a,b,c and d". These "parameters" are actually constant functions defined over s which is a distance.

These can be interpreted as the contents (heights) of some pollutants in the cross section of a river, whose length of interest is "s". The interactions of these pollutants gives rise to the mixing functions y(1) and y(2) which I'm trying to find.

thanks for the help!

Tags for this Thread

What are tags?

A tag is like a keyword or category label associated with each thread. Tags make it easier for you to find threads of interest.

Anyone can tag a thread. Tags are public and visible to everyone.

Contact us