Prev: simulation of refrigerator
Next: How to plot compass off origin?
From: Jan Simon on 4 Feb 2010 10:37
Dear Alexander!

> I believe I know how to fix 1. and 2. But concerning esp. 4: if it is not too comprehensive/does not involve much research, can you give me some more instructions what to do here? I took a look at the Vectorized section of the odeset help page, but actually didn't quite understand what to do and how to set that vectorized property. Can you give me some more instructions here?

The vectorized ODEFCN is an advantage, if you calculate the sensitivities of the result to the start values. This is not done in ODE45, but e.g. in ODE15S. The modifications are easy:
The ODE function must accept the input y and p as matrix and reply a matrix also. So the ode is evaluated once for the different columns of the input.

> Another thing: You recommended fast C integrator. However, I am much more fluent in FORTRAN as I recently did a project in this language. Do you happen to know a FORTRAN integrator which would be very fast?

"Very fast" is a relative statement for integrating. It is always possible to be very fast, but with a reduced precision and accuracy. Depending on the problem the result can be valid or nonsense! RK45 is a well method if the problem is not stiff, and there are lots of Fortran implementations in the net (I admit, I have no link, but Google has).
If you are doing the integration for scientific purposes, I'd suggest to include an analysis of the sensitivities and stiffness ever. A result is not a result, if the variation of the inputs by EPS effects a completely different solution! E.g. simulating the solar system for the next 10 years is valid if you consider the sun and the planets - but for the next 10 million years a "very fast" integration is meaningless if you do not consider the initial positions of the planets with and precision of cenitmeters, the moons and the interplanetary gas...

Good luck, Jan
From: Alexander Erlich on 5 Feb 2010 07:16
Thanks again Jan,

I did my best to implement 1., 2. and 4., but there are still some problems. I have pasted the new m-files below, but here are the essential changes:

In the DuffingODE file, I changed:
xdot=[x(2,:);-x(3,:).*x(2,:)+x(5,:).*x(1,:)-x(4,:).*x(1,:).^3+x(6,:).*cos(t);0;0;0;0];

In the Duffing_attractors file, I changed the ODE call to this:
[t,Xcurrent]=ode15s(@DuffingODE,interval.*omega_E, [x;y;C;alpha;beta;gamma]);

Your first two points were easy to correct: I expanded the initial conditions vector by the necessary parameters C, alpha, beta and gamma. Also, I removed omega_E from the initial conditions list and now integrate cos(t) over interval.*omega_E.

But concerning the vectorization, I think I would require several more instructions. I tried to keep as closely as possible to the example given in help odeset ---> Vectorized. However, my sort of vectorized version (which now uses ode15s) is considerably slower than the non-vectorized ode45. I suppose my vectorization is wrong here. Can you help me here?

Alexander

-------------------------DuffingODE.m------------------
function xdot = DuffingODE(t,x)
%DuffingODE contains set of two ODEs for Duffing oscillator.

xdot=[x(2,:);-x(3,:).*x(2,:)+x(5,:).*x(1,:)-x(4,:).*x(1,:).^3+x(6,:).*cos(t);0;0;0;0];

------------------------Duffing_attractors.m---------
function Duffing_attractors(interval, gamma, tol)
%DUFFING_attractors Analyzes the attractor basins of the Duffing oscillator.
%
% Example call: Duffing_attractors([0, 100],0,0.1);
%
% The corresponding ODEs are in the file DuffingDGL.m. They are taken from
% Argyris / Faust / Haase: Die Erforschung des Chaos. Vieweg, 1994, p. 669.
% The parameters (C, beta, alpha, gamma, omega_E) in the example call are
% chosen to be the same as in Argyris, p. 673.

%global C beta alpha gamma omega_E;
C=0.12; beta=0.5; alpha=0.5; omega_E=0.7;

clc; clf; clear All; % clear command, figure, all variables

figure(1);

h=plot(1,1); hold on;

for x=-4:0.1:4

if(ishandle(h)==0)
break;
end

for y=-2:0.1:2

if(ishandle(h)==0) % checks if figure has been closed. If so, ...

break; % ... simply leave the loop
end

[t,Xcurrent]=ode15s(@DuffingODE,interval.*omega_E, [x;y;C;alpha;beta;gamma]);
x1=Xcurrent(:,1);
% x2=Xcurrent(:,2); % x2 = dx1/dt will not be required.

% if long-term behaviour is attracted to 1 (i.e. [1,0] in the
% phase space portrait), draw a blue dot.

if(abs(x1(end)-1)<tol)
plot(x,y,'o','Color','Blue','linewidth',6)
drawnow; %update the plot
end

% if long-term behaviour is attracted to -1 (i.e. [-1,0] in the
% phase space portrait), draw a yellow dot.

if(abs(x1(end)+1)<tol)
plot(x,y,'o','Color','Yellow','linewidth',6)
drawnow; %update the plot
end
end
end
First  |  Prev  | 
Pages: 1 2
Prev: simulation of refrigerator
Next: How to plot compass off origin?