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.

6 个评论

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??
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
Jan 2017-8-27
编辑:Jan 2017-8-27
@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.
@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
@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.
@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.

请先登录,再进行评论。

回答(1 个)

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.

5 个评论

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.
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?
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'?
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".
Thank you so much for your help. I greatly appreciate it.

请先登录,再进行评论。

类别

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by