Help With ODEs Solver(ode45)

조회 수: 1 (최근 30일)
Ahmed Sulaiman Ahmed Al Rawahi
Ahmed Sulaiman Ahmed Al Rawahi 2021년 3월 12일
편집: David Goodmanson 2021년 3월 12일
function [] = call_func
clear all
close all
tspan = [ 0 10 ]
y0 = 0:0.01:1;
x0 = 0:0.01:1;
theta0 = 0:pi/5:2*pi;
l = 1;
options=odeset('RelTol',1e-8,'AbsTol',1e-8)
for k = 1:length(y0)
for s = 1:length(x0)
for w = 1:length(theta0)
[t,x] = ode45(@func,tspan,[x0(s) y0(k) theta0(w)],options);
%legend('b = blue', 'r = red')
if (max(mod(x(:,1),l))<0.5) && (min(mod(x(:,2),l))>0.5)
c='r.';
else
c='b.';
end
plot(mod(x(end,1),l),mod(x(end,2),l),c);
hold on
end
end
end
end
function dydt = func(t,x)
Vs = 0.095;% true value
alpha =1; % true value
dydt = [sin(2*pi*x(1))*cos(2*pi*x(2))+Vs*cos(x(3))
-cos(2*pi*x(1))*sin(2*pi*x(2))+Vs*sin(x(3))
2*pi*(-alpha*cos(2*pi*x(1))*cos(2*pi*x(2))*sin(x(3))+sin(2*pi*x(1))*sin(2*pi*x(2)))];
end
I use this code to analyze a fluid system. The problem is it take sometimes up to 20 mintues to get results. Does anyone know if there is a way to solve and get the plot a quicker way?

채택된 답변

David Goodmanson
David Goodmanson 2021년 3월 12일
편집: David Goodmanson 2021년 3월 12일
Hi Ahmed,
repeating the plot command every time in the nested for loops is not helping, For a smaller case (I didn't want to wait for 20 minutes), the following code plots once and cuts the time down by about a factor of 5:
clear all
close all
tspan = [ 0 10 ]
y0 = 0:0.01:.1;
x0 = 0:0.01:.1;
theta0 = 0:pi/5:2*pi;
l = 1;
options=odeset('RelTol',1e-8,'AbsTol',1e-8);
Xr = zeros(length(x0),length(y0),length(theta0));
Yr = zeros(length(x0),length(y0),length(theta0));
Xb = zeros(length(x0),length(y0),length(theta0));
Yb = zeros(length(x0),length(y0),length(theta0));
%tic
for k = 1:length(y0)
for s = 1:length(x0)
for w = 1:length(theta0)
[t,x] = ode45(@func,tspan,[x0(s) y0(k) theta0(w)],options);
%legend('b = blue', 'r = red')
if (max(mod(x(:,1),l))<0.5) && (min(mod(x(:,2),l))>0.5)
Xr(s,k,w) = mod(x(end,1),l);
Yr(s,k,w) = mod(x(end,2),l);
else
Xb(s,k,w) = mod(x(end,1),l);
Yb(s,k,w) = mod(x(end,2),l);
end
% fig(1)
% plot(mod(x(end,1),l),mod(x(end,2),l),c);
% hold on
end
end
end
%hold off
fig(2)
plot(Xb(:),Yb(:),'b.')
%toc
  댓글 수: 2
Ahmed Sulaiman Ahmed Al Rawahi
Ahmed Sulaiman Ahmed Al Rawahi 2021년 3월 12일
Hello David.
I appriecate your help! That is much faster and gave me a better results.
If you have time, would you mind walking me through the changes you made and brefly explain what you did ?
Thank you again.
David Goodmanson
David Goodmanson 2021년 3월 12일
편집: David Goodmanson 2021년 3월 12일
Hi Ahmed
You have three nested for loops with indices k,s,w, so for each occurrance in the for loops, you can save the X and Y data as 3 dimensional arrays with indices k,s,w. There are separate arrays for red and blue (does red occur?)
After all the calculations are done, for plotting purposes the 3d arrays won't work. The expression Xb(:) reads out, in a particular order, the 3d array into one long column of data. Yb(:) also reads out into a column, of course in the same order, so the Xb and Yb values match up. With those two columns you can make an xy plot that matches the plot-every-time procedure.
Preallocation of the 3d arrays occurs up toward the top.

댓글을 달려면 로그인하십시오.

추가 답변 (0개)

카테고리

Help CenterFile Exchange에서 Loops and Conditional Statements에 대해 자세히 알아보기

Community Treasure Hunt

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

Start Hunting!

Translated by