开发者

Finding Intersection Points between 3rd Order ODE and a line?

开发者 https://www.devze.com 2023-02-08 22:44 出处:网络
How can I find the intersection poi开发者_C百科nts between the solution of a 3rd order ODE and a line y=x?

How can I find the intersection poi开发者_C百科nts between the solution of a 3rd order ODE and a line y=x?

My ODE's code is

sol=dsolve('D3y-4*D2y+Dy+2*y=0,y(0)=-4,Dy(0)=-6,D2y(0)=-4')
x=0:2
y=subs(sol,'t',x)
plot(x,y)


To find the numerical values of the intersections:

Create an anonymous function that is zero at the intersections between your solution and the line you seek:

sol = dsolve('D3y-4*D2y+Dy+2*y=0,y(0)=-4,Dy(0)=-6,D2y(0)=-4');
my_func = @(x) subs(sol,'t',x) - x; % Your solution - x is equal to zero at the %intersections

then, if you want to find the values graphically:

plot(-5:.01:5,my_func(-5:.01:5))

or numerically, by optimization routine:

x = fzero(my_func,0); % I find x = -.6847

which will look for a zero of the function near to 0. It will not find all the zeros so you will need to start the fzero function near values where you expect the intersections to be.

Hope this helps,

Andrew

Edits: a how to of the bisection method.

Once we have our anonymous one dimensional function "my_func" if you don't want to use an optimization method to solve the equation, but happen to know a range [range_{min} range_{max}]in which my_func = 0 then the following algorithm will find the zero of the function for you, provided you're working with a continuous function:

range_min = 0; % say our range is [0 2]
range_max = 2;

error_tolerance = .0001; % will find the answer to within .0001

while (range_max - range_min < error_tolerance)
range_temp = (range_max + range_min)/2;
if ((my_func(range_temp) <0 & my_func(range_max)>0) | (my_func(range_temp) >0 & my_func(range_max)<0))
range_min = range_temp;

else if ((my_func(range_min)<0 & my_func(range_temp)>0) | (my_func(range_min)>0 & my_func(range_temp)<0))
    range_max = range_temp;

else if (my_func(range_temp == 0)
range_min = range_temp;
range_max = range_temp;

    end

end

t_intersection = (range_min + range_max)/2;

So some explanation: if the function is continuous, and yours is, then if it intersects y=t at t_intersection then at our modified function my_func(t) = sol(t)-t will have a zero at t_intersection. Since my_func is continuous then we can find a zero of the function provided we know two values of the function, one that is greater than zero and one that is less than zero.

So we start with these known points and define a range [range_min range_max] where either my_func(range_min)<0 and my_func(range_max)>0 or vise-versa. Then we cut this range in half, by creating the midpoint range_temp = mean(range_min and _max). We now create a new range [range_temp range_max] or [range_min range_temp] so that we preserve a sign change of my_func over the range. We repeat this process until we reach a satisfactory accuracy.

One caveat though, this will only find one zero within the initial range you provide. That is a fundamental frustration with most zero finding methods, and more generally the field of optimization of which zero finding can be considered a special case.

  • If you need more on the bisection method: http://en.wikipedia.org/wiki/Bisection_method
  • If you need more on zero finding methods, I would look in a university library for a book on numerical analysis.
  • If you need more on anonymous functions in Matlab (the " = @(t)" thing I was using), then check here: http://www.mathworks.com/help/techdoc/matlab_prog/f4-70115.html

I think that covers it, good luck.

--Andrew

0

精彩评论

暂无评论...
验证码 换一张
取 消