Parameter calculation by using loops in ODE45
1 次查看(过去 30 天)
显示 更早的评论
For any initial values(Ofcourse with limits) I used it either excutes on the 4th iteration or it would be beyond the MAX iteration.
I am afraid something is wrong and I need a help.
zspan=[0,400];
v0mat = [1 0.05 1];
tol = 10^-4; % Treshold
MAX = 1000;
v2 = 0.4;
interval = [0.01 0.2]; % alpha interval
a = interval(1);
b = interval(2);
alpha = (a+b)/2;
zsol = {};
v1sol = {};
v2sol = {};
v3sol = {};
for k=1:size(v0mat,1)
v0=v0mat(k,:);
[z,v]=ode45(@(z,v)rhs(z,v,alpha),zspan,v0);
zsol{k}=z;
v1sol{k}=v(:,1);
v2sol{k}=v(:,2);
v3sol{k}=v(:,3);
end
iter = 1;
while(abs((v(:, 2)) - v2) > tol)
alpha= (a + b)/2;
[z,v]=ode45(@(z,v)rhs(z,v,alpha),interval,v0);
if(abs(v(:,2))-v2 > 0)
b = alpha;
else
a = alpha;
end
iter = iter + 1;
if(iter > MAX)
return;
end
end
for k=1:size(v0mat,1)
figure(1)
plot(v2sol{k},zsol{k},'g')
xlabel('Velocity,w')
ylabel('Height, z')
grid on
end
function parameters=rhs(z,v,alpha)
Nsqr = 0.001;
db= 2*alpha-(v(1).*v(3))./(2*v(2).^2);
dw= (v(3)./v(2))-(2*alpha*v(2)./v(1));
dgmark= -Nsqr-(2*alpha*v(3)./v(1));
parameters=[db;dw;dgmark];
end
16 个评论
Walter Roberson
2020-1-17
编辑:Walter Roberson
2020-1-17
In your context, timepoint is z value. Your ode45 call will process over z in zspan, producing one row for each of a number of z values. (abs(v(:,2))-v2 > 0) is a vector result that might be true for some of the z entries but false for other times. You need to decide which z the condition must hold for.
回答(0 个)
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Loops and Conditional Statements 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!