Question about using ODE45

I have this equation: dy/dx= 1/(a*x^-3+b*x^-5)

a=(3/4)t

b=(1/2)t^2

t=[0.5 20]

x=[2 5]

The question is using ODE45 to solve that equation for every value of t

This my function file:

function [dydx] = Myfun( x,y,t )
dydx=1/(((3/4.*t).*x^(-3))+((1/2).*t^2).*x^(-5))
end

My script file:

t=linspace (0.5,20,10);
options=optimset('Display','Off');
for in = 1:length(t)   
  T1=t(in)
  [x,y,x1,y1,i1]=ode45(@(x,y)Myfunc(x,y,t(in)),[2 5],0,options);
  xn=x1(:,1)
  u(in)=x1(:,1);  
end

For some reason, my script doesn't run and it shows error. Please kindly help me to check what the problem is.

Thanks in advance

Answer 1

[Star Strider]

You need to vectorise ‘Myfunc’. You also need to check the initial condition, since a zero initial condition results in a uniformly zero integrated result.

Try this variation on your code:

Myfunc = @(t,x) 1./(((3/4.*t).*x.^(-3))+((1/2).*t.^2).*x.^(-5));

t = linspace (0.5,20,10);
options=optimset('Display','Off');
[T,X] = ode45(Myfunc, t, 1E-3, options);

figure(1)
plot(T, X)
grid on

[Nicia Nanami]

Thank you for your help. I appreciate it!!

However the question asks me to find y for every of t. Thus I suppose the answer must like t=0.5, y=??. I have changed my function file a bit:

[ dydx ] = Myfun( X,y,t )

dydx=1/(((3/4.*t).*X(1)^(-3))+((1/2).*t^2).*X(1)^(-5));

And this is my scipt:

t=linspace (0.5,20,10);

options=optimset('Display','Off');

for in = 1:length(t)  

    T1=t(in)

    [x,y,x1,y1,i1]=ode45(@(X,y)Myfun(X,y,t(in)),[2 5],0,options);

    yn=y1(:,1)

end

But I still couldn't make it run properly.

[Nicia Nanami]

It states "Output argument "varargout{4}" (and maybe others) not assigned during call to "ode45".

Thank you

[Walter Roberson]

The TE, YE, and IE outputs (ode45 outputs #3, #4, and #5) are only defined if you have set an event function in options.

[Star Strider]

@Nicia Nanami — My pleasure. Walter’s Comment completely addresses your question about the varargout error. Use my code and it will not appear. My code also calculates ‘y for every t’. You absolutely do not need the loop!

@Walter — Thank you!

[Nicia Nanami]

@Walter: thank you for your explanation.

@Star Strider: thank you so much for your help. I have tried to use your code but the outputs are t and x not t and y

[Star Strider]

@Nicia Nanami — My pleasure.

The outputs in my ode45 call are ‘T’ and ‘X’!

Please read — and understand — the code I posted!

[Nicia Nanami]

I got it. Thank you so much for your help

[Nicia Nanami]

I have one short question. In this question, the range of x is from 2 to 5. So how can I apply this condition into your code?? I'm so sorry if any convenience that I may cause to you.

[Star Strider]

My pleasure.

You can supply the initial condition (third argument to ode45) as 2.

Myfunc = @(t,x) 1./(((3/4.*t).*x.^(-3))+((1/2).*t.^2).*x.^(-5));

t = [0.5  20];
options=optimset('Display','Off');
ic = 2;
[T,X] = ode45(Myfunc, t, ic, options);
T(X >= 5) = [];
X(X >= 5) = [];

figure(1)
plot(T, X)
grid on

You could probably set an event function to stop the integration when X>=5. Setting all those values to the empty array ‘[]’ is just easier.

No apology necessary!

[Nicia Nanami]

Thank you for your help, I appreciate it.

[Star Strider]

As always, my pleasure!

Answer 2

[Roger Stafford]

(Corrected)

Using ‘ode45’ is a very bad idea for this problem. Since dy/dx does not depend on y, this is a simple problem in integration for which you can use Matlab’s ‘int’ function in the Symbolic Toolbox. You can get a solution from ‘int’ in terms of 'a' and ‘b’, taking into account your y value of zero for x equal to 2. All you have to do then is substitute in the respective values 3/4*t and 1/2*t^2 and you have a general expression for y in terms of x and t.

[Nicia Nanami]

Thank you for your respond but the question asks to use ode45 to solve this problem so I'm allowed to use other technique to solve it.

[Jan]

@Nicia: I assume this is a homework of a lesson for numerical mathematics. Then Roger's argument is essential (+1), because actually the instructions should teach you how to use numerics properly. Using ODE45 to integrate this is a poor approach. I have seen too many works of students and PhD theses, which used the wrong numerical method with the rationale, that they "produce a result". It is like using a hammer to push a screw into to wood.

But this is a side note only. If your professor asks for using ODE45, this is his mistake, but you should use it for this homework.

Answer 3

[John BG]

Hi Nicia

this is John BG <mailto:[e-mail] [e-mail]>

1.

t is just a parameter, not the reference vector of the differential equation.

the output of ode45() is not

[T,X]=ode45(..)

but

[X,Y]=ode45(..)

.

2.

There are other more efficient ways to use ode45, but a way to have it running is

clear all;clc;close all;

x=[2:.01:5];
t=[-100 .5 20 1e3];

options=optimset('Display','off');

x_span=[2 5];

figure(1); hold all

for k=1:1:numel(t)
    t1=t(k);
    f1=@(t1,x) 1/((.75*t1)*x^-3+(.5*t1^2)*x^-5);
    [X,Y]=ode45(f1,x_span,0.7,options)  
    hold all; plot(X,Y)
end

grid on

.

3.

Perhaps a 'better picture' of the beahviour is obtained when expanding the range

x=[-20:.01:50];
x_span=[-20 50];

figure(2); hold all

for k=1:1:numel(t)
    t1=t(k);
    f1=@(t1,x) 1/((.75*t1)*x^-3+(.5*t1^2)*x^-5);
    [X,Y]=ode45(f1,x_span,0.7,options)  
    hold all; plot(X,Y)
end

grid on

.

4.

Also, integration required, because of the dydx, than then would be

clear all;clc;close all;

x=[2:.01:5];
t=[-100 .5 20 1e3];

options=optimset('Display','off');

x_span=[2 5];

figure(1); hold all

for k=1:1:numel(t)
    t1=t(k);
    f1=@(t1,x) 1/((.75*t1)*x^-3+(.5*t1^2)*x^-5);
    [X,Y]=ode45(f1,x_span,0.7,options)  
    hold all; plot(X,cumsmum(Y))
end

grid on

5.

and integrating the wider range

x=[-20:.01:50];
x_span=[-20 50];

figure(2); hold all

for k=1:1:numel(t)
    t1=t(k);
    f1=@(t1,x) 1/((.75*t1)*x^-3+(.5*t1^2)*x^-5);
    [X,Y]=ode45(f1,x_span,0.7,options)  
    hold all; plot(X,cumsum(Y))
end

grid on

.

.

if you find this answer useful would you please be so kind to consider marking my answer as Accepted Answer?

To any other reader, if you find this answer useful please consider clicking on the thumbs-up vote link

thanks in advance

John BG

<mailto:[e-mail] [e-mail]>

[Nicia Nanami]

Thank you so much, your work is really clear and helpful. But I have one question. If I want to plot t vs y, how can I do it??

[Jan]

@John BG: This is not useful:

for k=1:1:numel(t)
  t1=t(k);
  f1=@(t1,x) 1/((.75*t1)*x^-3+(.5*t1^2)*x^-5);
  ...
end

Your definition of f1 does not depend on the value of t1=t(k). Therefore your loop for k=1:1:numel(t) is meaningless, because the body does not depend on the loop counter k. Perhaps Nicia has overseen this.

Your points 3 and 5: It is not meaningful to let ODE45 integrate over the discontinuity near x=0. ODE45 cannot handle discontinuities properly and the replied trajectory is junk. The original question defines "x=[2 5]", such that the behavior near to 0 was not part of of the question.

[John BG]

Hi again Nicia

ok, please confirm; you need t(y), not y(t) right?

would it be possible to broadly know what are you modelling, a rupture?

[John BG]

Jan Simon, would you please abstain from adding any comment to my question to Nicia t(y) or y(t)? at least until Nicia says something?

Nicia found my answer helful, what was the point of you saying not useful? it's not your question aber Nicia's.

If Nicia wishing to do so, please let Nicia follow up right after my question.

And Jan, you didn't span the range to see the step, I did show the step. Spotting when it happens is probably the objective of the question.

Again, To Nicia:

Hi again Nicia

do you need t(y) or y(t)?