Parametrizing ode45

Hi everybody,

I am new to Matlab and would request some help with numerically solving an ODE. The equations themselves are very simple as given below:

u1'=-u1*u2*cos(theta) u2'=-u1^2*cos(theta)

Where u1 is u1(z) and u2 is u2(z) and the derivative is w.r.t z. I would like to generate plots for z=[0,5] and theta=[-2pi,2pi]. I tried the simple case where theta is fixed to pi and obtained a nice plot for u1,u2 vs z.

function dydz=amstrong(z,y)
theta=pi;
dydz=zeros(size(y));
u1=y(1);
u2=y(2);
dydz(1)=-u1*u2*cos(theta);
dydz(2)=u1^2*cos(theta);
end

I would like to do this for a range of theta values. I am not sure how to proceed.

The script to get results is as shown below:

    %%%%%%%
    tic
    %%%%%%%
    clear % Clear memory

    z0=(0:0.0001:5);                % Vector that goes from 0 to 5 
    y0=[1,0];                       % IC u1=1 and u2=0    

    %Calling ode45
    %
    [z,y]=ode45('amstrong',z0,y0);

    %Plotting the result.The plot is of u^2 v/s z.
    plot(z,y.^2) 

    %%%%%%%%
    toc
    %%%%%%%%

My objective to get a series of graphs on the same plot as theta is varied over a large range. I did take a look at the possibility of a function calling another function, but it was very confusing.

Any help will be much appreciated.

Thank you.

PS# Running Matlab2010a x64 on Windows7.

Answer 1

[Jan]

See this thread, where Jiro explains using anoymous function for parameterization:

Answers: What is the benefit of anonymous functions

Answer 2

[Abhishek Murthy]

1. declare theta as global variable in both the functions.

2. In the "script to get results", vary theta in a loop. Inside the loop call ode45. In every loop iteration, armstrong will be called upon to provide the system definition. It would read the current value of theta as its globally visible.

3. process [z,y] vectors in each loop into some data structure (maybe a 3d matrix or a linked list) where the third dimension indexes various values of theta

I hope this helps.

Answer 3

[Antillar]

Abhishek: Thank you for the algorithm. I will get back to you after giving it a go.

Jan: Thanks for the link. I'm still digesting Jiro's post. Any restrictions on the type of object I can pass as a parameter (i.e scalar,vector,matrix?)

Thank you guys.

[Jan]

You can pass anything as parameter to anonymous functions, even further anonymous functions.

Answer 4

[Antillar]

Ok guys this seems to work, but is exceedingly slow. :(

% Defining the function SHG 
  function [du_dz]=SHG(z,u,theta)
  du_dz=zeros(size(u));
  u1=u(1);
  u2=u(2);
  du_dz(1)=-u1*u2*cos(theta);
  du_dz(2)=u1^2*cos(theta);
  end

I use the function as follows:

% Solution to:
% du_dz(1)=-u(1)*u(2)*cos(theta)
% du_dz(2)=-(u(1)^2)*cos(theta)

% set an error
options=odeset('RelTol',1e-12); 

% Set the span of z and theta and get the number of columns.
z=(0:1:5);
[row1,col1]=size(z);
theta=(-2*pi:0.1:2*pi);
[row2,col2]=size(theta);

% Initial conditions
u0=[1;0];

% Initialize Matrix x. This matrix contains values of u1 and u2 at each theta. It is an n:2:m dimensional matrix, with n as the columns of z and m the columns of theta as the third dimension.
x=zeros(col1,2);

% Call the solver by passing thetas as a parameter.
for j=1:col2,
[z,u] = ode45(@SHG,z,u0,options,theta(j));
x=cat(3,x,u); %Expand the Matrix in the 3rd dimension
end

%Plot the results of u1^2 and u2^2 v/s theta for a user supplied z.

clear j;
z_value=input('Enter a z value: '); %User supplies a value of z
index=find(z==z_value); %Get the index of the array element
for j=1:col2,
%figure 
hold on
plot(theta,x(index,1,j).^2)
plot(theta,x(index,2,j).^2)
   legend('Fundamental','Harmonic');ylabel('Intensity');xlabel('phase')

end

Did is I mess something up because it is taking a ridiculously long time to plot the graphs. I had to kill the program.

Any help will be appreciated!

Thank you.