Solve ODES with multiple initial conditions

global tau;
global T0; To=120+298
global V_R; V_R=2*1e-5
function dy=ODES(z,y)
%ODES
dy(1)=V_R*(-2*r(1)-2*r(2));
dy(2)=V_R*(-0.5*r(1));
dy(3)=V_R*(r(1)-r(3));
dy(4)=-V_R*(DELH1*R(1))/y(1)/CPA
dy=dy';
dy=dy(:);
end
CT0=100000/8.3145/T0;  
v0_0=V_R/tau;         
FT0=v0_0*CT0;       
R0=[8 1 0]/(8+1);
C0=R0*CT0;
F0=v0_0*C0;

tspan=[0,1];
y0=[F0,T0];
[z,y]=ode15s(@ODES,tspan,y0);

I want to sovle the ODES with differen tau from 0.1 to 20 and than plot the tau versus y(2) at the z=1. Does someone know how to do it.

Appreciate your help

[Walter Roberson]

No, we do not know how to do it.

Your y0 depends upon T0 which is not defined.

Your y0 depends upon F0 and T0. F0 is length 3. T0 is unknown length. The number of outputs from the ODE must match the number of entries in y0, so with your y0 being 3 plus some unknown length, your ODES must return 3 plus that unknown length, but instead it only returns length 3.

Your ODES depend upon V_R and r(1) and r(2) and r(6) but V_R and r are not defined.

Your ODES do not depend upon the current boundary conditions in y.

There is a R0 in the initialization but that creates a vector of length 3 and so does not match up to r(6) being needed.

[Zhihong Lin]

Hello Walter,

I am sorry that these code makes you confused. i just randomly delete lots of codes because the original code is kind of long to post. could you help me how to vary the initial conditions so that i could plot like what i said before. thank you very much.

[Walter Roberson]

Where do you define r and R that you use in the function? And DELH1 and CPA? Where do you initialize T0? You initialize To which is a different variable.

Answer 1

[Athul Prakash]

Hey Zhihong Lin,

As has been pointed out by others, your initialization code is a bit confusing and buggy. I understand that you posted only a shortened version here. Hence, I am assuming that, given any value of tau, you have the code to generate all the required parameters for the Diff Equations.

Three things I think you should consider:

1) For-Loop:

Have you tried a simple for loop to iterate through your range of tau; each time calling ode15s to solve the problem for a different initial condition. Here’s what I mean, through pseudo-code:

tau_vector = 0.1:step:20;   %choose your step value
for tau_i = tau_vector:
    % Here, Initialize all your variables that depend on tau_i.
    [z, y] = ode15s(@ODES, tspan, y0);
    % Here, add y(2) into a vector which would hold the values of all iterations.
end
% now you can plot the record of y(2) against tau_vector.

Ofcourse, if you have too many values in tau_vector, this approach might take a long time. In that case, consider the following methods to speed it up.

2. Use parfor

parfor, from the Parallel Computing Toolbox, lets you run multiple iterations of a for loop in parallel (on different CPU cores). Refer to the documentation below:

https://www.mathworks.com/help/parallel-computing/parfor.html

The iterations of the loop are independent of each other, hence parfor is suitable here, if required.

3. Vectorize all your variables into 1 ODE System

If you have N values of tau, treat your differential equations as a 4xN variable problem.

Instead of 4 variables, y(0) to y(4), use y(0) to y(4*N). Vectorize all the variables required and use ode45 to solve all equations simultaneously, as if they were a single system of 4*N equations.

Depending on your exact code, this might take an even longer time though. But I think it’s worth trying out.

Good Luck!