4 views (last 30 days)
Nicia Nanami on 4 Oct 2017
Edited: Jan on 10 Oct 2017
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.

Star Strider on 4 Oct 2017
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
Star Strider on 6 Oct 2017
As always, my pleasure!

Roger Stafford on 4 Oct 2017
Edited: Roger Stafford on 4 Oct 2017
(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.
##### 2 CommentsShowHide 1 older comment
Jan on 6 Oct 2017
Edited: Jan on 6 Oct 2017
@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.

John BG on 5 Oct 2017
Hi Nicia
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
.
.
John BG
John BG on 7 Oct 2017
Edited: John BG on 7 Oct 2017
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)?