Please help me with the code

I have an equation:

 dy/dx = 8*y^k - 3*y

The given numbers:

k=[0.1:0.01:0.2]
y0=0.6
x=[0 10]
y1=2

How can I use ODE to solve the equation and use that solution to find x1 (x1 is at y1) for every value of k. CANNOT use arrayfun or find.

[KALYAN ACHARJYA]

here x=1/(8y^(k-1)-3), Now put the K values 0.1000 0.1100 0.1200 0.1300 0.1400 0.1500 0.1600 0.1700 0.1800 0.1900 0.2000... After that question is not clear. y0=0.6??

[Image Analyst]

Why can't you use arrayfun or find()? Usually these kinds of weird restrictions are only for homework problems. But you didn't tag it as homework. Is it? http://www.mathworks.com/matlabcentral/answers/8626-how-do-i-get-help-on-homework-questions-on-matlab-answers

[Jan]

@Nicia: "x1 is at y1" is not clear to me. Actually all you need is a loop over the elements of k and a call of ODE45. But the "y1=2" is confusing. Is this an initial or boundary value problem?

If this is a homework, I will not post the complete solution to give you a chance to submit your solution.

[Nicia Nanami]

@KALAN, @Jan Simon: Thank you for replying my post. This question asks to find x1, when y is equal 2, it's called y1. At that time, x is equal x1. I understand how to use ODE45 to find the solution of the equation. This is what I do. I create an m.file

[ dydx ] = func( x,y)

k=0.1

dydx=(8.*(y.^k))-(3.*y);

[x,y]=ode45(@func,[0:0.5:10],0.6)

The output will be y over x. However, I don't know how to find x1 without using find or arrayfun function. Do you know any functions that I can use in this situation? And since k is a range that means I have to create many m.file for each k value. Is there any better ways for me, instead? I just want to be given a hint not a complete solution. Thank you in advance

[Nicia Nanami]

@Image Analyst: yes this is my homework. I'm so sorry for not tagging it as homework since this is my first time using this website. Thank you for reminding me.

[Jan]

@Nicia: This looks confused:

[ dydx ] = func( x,y)
k=0.1
dydx = (8.*(y.^k))-(3.*y);
[x,y]=ode45(@func,[0:0.5:10],0.6)

Better:

k     = 0.1
dy    = @(x,y) 8 * (y.^k) - 3 * y;
[x,y] = ode45(@dy, 0:0.5:10, 0.6)

Remark: Unnecessary square brackets should be avoided (see Why not use square brackets). I prefer the use of optional round parentheses, if they improve the clarity of the code. A multiplication has a higher precedence than the subtraction, so here I would omit the parentheses, but e.g. -2^k could be confused with (-2)^k, so I write -(2^k). It's a question of taste.

Answer 1

[Jan]

If this is an initial value problem, all you need is a loop over the elements of k:

kList = 0.1:0.01:0.2;

Now pick a specific k and:

k      = kList(1);
dy     = @(t, y) 8*y^k - 3*y;
[x, y] = ode45(dy, [0, 10], 0.6);

Is the last element of y the searched value? Then all you need to do is to embed this in a loop and collect the results.

[Nicia Nanami]

Thank you for helping me. The question asks to find x1( x is equal x1 when y reaches 2 called y1). I have found the solution of the equation. The output is y over x. But I'm having difficult time to find x1 without using find or arrayfunc function. Please guide me what I should do next. Thank you in advance.

[Jan]

Does this mean that you want to stop the integration when y=2 is reached? Then an event function will help: See doc odeset. fzero might be useful also. I have no idea, why arrayfun could be used here.

Do not create different M-files for the elements of k, but a loop:

kList = 0.1:0.01:0.2;
for ik = 1:length(kList)
  k = kList(ik);
  dy     = @(t, y) 8*y^k - 3*y;
  ...
end

Is this a homework for numerics and you have heard something about boundary value problems with free end time?

[Nicia Nanami]

I haven't heard anything about boundary value problems with free and time. What I intend to do is like drawing out x1 from the solution got by using ode45 to solve the equation. And x1 satisfies a condition that y=2. I don't if any function can help me with that. I have been searching around but I couldn't find it. That's why I come here and ask for help. However, there maybe another way to solve this problem better than mine.

In your code: dy=@(t,y).. what is 't', does you mean 'x'?

[Jan]

Yes, the "t" could be a "x" also. But as long as both do not appear in the calculations, this does not matter. < must be the 2nd input, that's all.

As said, you can get the x value, at which y=2 by using an event function for the integrator. This is defined by odeset. This function can stop the integration when y=2 is reached and then you get the corresponding x value as output also:

opt    = odeset('Events', @myEventsFcn);
[x, y] = ode45(dy, [0, 10], 0.6, opt);

function [value, isterminal, direction] = myEventsFcn(x, y)
...

Now check for y==2 and set the isterminal flag to stop the integration.

See: doc ode45 and search for "ODE Event Location".

[Nicia Nanami]

Thank you so much for your help. I greatly appreciate it.