Determine where lines intersect

1 回表示 (過去 30 日間)
Andreas
Andreas 2014 年 9 月 29 日
コメント済み: Jan 2014 年 10 月 5 日
Hey! Im wondering how i could Determinate the intersect of these two lines. I got the following code.
z0=10000;
dz0=0;
span=[0 250];
[T,Y] = ode45(@function_A,span,[z0 dz0]);
x=linspace(0,300,2000);
plot(T,Y(:,1),'b', x,0,'g');
and this function
function dz = function_A(t,z)
dz = zeros(2,1);
g=9.81;
m=250;
lambda=7.5*10^3;
k0=(m*g)/((250/3.6)^2);
k=@(z)k0.*exp(-z./lambda);
dz(1)=z(2);
dz(2)=-g-(k(z(1)).*abs(z(2))*z(2))/m;
end
Anyone who can help me?

回答 (1 件)

Mischa Kim
Mischa Kim 2014 年 9 月 29 日
編集済み: Mischa Kim 2014 年 9 月 29 日
Hello Andreas, use ode events. This examples shows how to detect events with the bouncing ball problem. Check out the following:
function my_ode()
z0 = 10000;
dz0 = 0;
span=[0 250];
options = odeset('Events',@events);
[T,Y,Te,Ye,Ie] = ode45(@function_A,span,[z0 dz0],options);
x = linspace(0,300,2000);
plot(T,Y(:,1),'b', x,0,'g');
end
function dz = function_A(t,z)
dz = zeros(2,1);
g = 9.81;
m = 250;
lambda = 7.5*10^3;
k0 = (m*g)/((250/3.6)^2);
k = @(z)k0.*exp(-z./lambda);
dz(1) = z(2);
dz(2) = -g-(k(z(1)).*abs(z(2))*z(2))/m;
end
function [value,isterminal,direction] = events(t,y)
% Locate the time when height passes through zero in a
% decreasing direction and stop integration.
value = y(1); % Detect height = 0
isterminal = 1; % Stop the integration
direction = 0; % Negative direction only
end
  1 件のコメント
Jan
Jan 2014 年 10 月 5 日
The function to be integrated contains an abs() . This is a non-smooth function and inconsequence the solution suffers from the problems described here. The step size control of ODE45 must fail at the changes of the sign. To handle this, the step size is reduced until the discontinuity is covered by the rounding errors.
Better use an event function and restart the integration at the discontinuity.

サインインしてコメントする。

カテゴリ

Help Center および File ExchangeOrdinary Differential Equations についてさらに検索

Community Treasure Hunt

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

Start Hunting!

Translated by