How to run a loop in linear differential equations in symbolic atmosphere?

Hi,

I have 6 symbolic functions i.e. y1 - y3 and c1 - c3, which I have used for defining linear ODE's.

The linear OD equations are as follows:

dy1/dt = (-c1*y1)+(c2*y2);

dy2/dt = (c3*y1)+(c2*y2)+(c3*y3);

dy3/dt = (c1*y1)+(c5*y2)+(c3*y3).

Now I want to define c's as constant values between two limits, i.e. c1=0.1 to 0.9; c2= 1 to 2.2 and c3=10 to 100.

Now what I want, to see effect of these equations with different constant values (c's), i.e. in loop form.

That is, I want to do two things:

1. Change c1 (from 0.1 to 0.9) and keep c2=1 and c3=10;

2. Change c1 (from 0.1 to 0.9) as well as c2 (from 1 to 2.2) and c3 (from 10 to 100).

And see values of y1, y2 and y3 in workspace or graph.

-----------------------------------------------------

I tried many times but all the time I get new type of error!! This time it is "In an assignment A(I) = B, the number of elements in B and I must be the same."!!

Please help me out. Thanks in advance.

Answer 1

[Star Strider]

The Symbolic Toolbox is not the best way to solve your problem. Use the numerical ODE solvers instead.

This will work, although you didn’t say anything about c5, so I assigned it a value of 1 to test the loops. I tested these to be sure they’d work, but I didn’t take the time to let the second loop complete. Be sure that the ‘tspan’ vector has the same length for all runs in a loop.

fun1 = @(t,y,c1,c2,c3,c5) [((-c1*y(1))+(c2*y(2)));  ((c3*y(1))+(c2*y(2))+(c3*y(3)));  ((c1*y(1))+(c5*y(2))+(c3*y(3)))];

% Loop #1:
T0 = clock;
tspan = linspace(0, 1, 25);
c1 = linspace(0, 0.9, 10);
c2 = 1;
c3 = 10;
c5 = 1;
for k1 = 1:length(c1)
    [t,y] = ode45(@(t,y) fun1(t,y,c1(k1),c2,c3,c5), tspan, ones(1,3)*eps);
    ty1(k1,:,:) = [t  y];
end
T1 = clock;
% Loop #2:
tspan = linspace(0, 1, 25);
c1 = linspace( 0, 0.9, 10);
c2 = linspace( 1, 2.2, 10);
c3 = linspace(10, 100, 10);
c5 = 1;
for k1 = 1:length(c1)
    for k2 = 1:length(c2)
        for k3 = 1: length(c3)
            [t,y] = ode45(@(t,y) fun1(t,y,c1(k1),c2(k2),c3(k3),c5), tspan, ones(1,3)*eps);
            ty2(k1,k2,k3,:,:) = [t  y];
        end
    end
end
T2 = clock;
ts = datestr(now, 'yyyymmdd-HHMMSS');
save(['LearnerODEmtxs_' ts '.mat'], 'ty1', 'ty2', 'T1', 'T2', 'T3')

This looks like it’s going to take a while, so I included a couple timers and a ‘save’ statement to store the results so you won’t have to run them again unless you change the conditions. The filename is unique for each run, and incorporates the date and time the file was created.

[Learner]

many many thanks for your help.

The loops is running. But there is one problem: It is not saving the results as you have mentioned. It is showing this error:

----------------------

Error using save Unable to write file LearnerODEmtxs_20140510-201923.mat: permission denied.

Error in Untitled (line 30) save(['LearnerODEmtxs_' ts '.mat'], 'ty1', 'ty2', 'T1', 'T2')

----------------------

ANOTHER thing I want to ask, can I see the effect of changing y's on graph (for both loops)? i.e. how changing c's changing the values of y's for defined tspan, on graph?

Thanks in advance.

[Star Strider]

My pleasure!

The save worked for me, but I just discovered that I copied an earlier version of the save line and didn’t catch the error. Replace it with:

save(['LearnerODEmtxs_' ts '.mat'], 'ty1', 'ty2', 'T0', 'T1', 'T2')

If that still gives you problems, then it has to do with your ability to save files on the computer you are using. (You can also delete or comment-out the ‘save’ line and the one before it if you don’t want to save the data.)

It may be difficult to figure out how to plot your results, especially for ty2. Here is a way to plot y(1) from ty1 as a function of c1 and time to get you started:

% Plotting ‘y(1)’ from ‘ty1’
[X,Y] = meshgrid(tspan,c1);
Cplot = squeeze(ty1(:,:,2:4));
figure(1)
surf(X, Y, squeeze(ty1(:,:,2)))
xlabel('Time')
ylabel('C_1')
zlabel('Y_1')
title('Y_1 as a function of C_1 and Time')
grid on

It seems c1 doesn’t have a noticeable effect. Note that y(2) is ty1(:,:,3) and y(3) is ty1(:,:,4). When I did some exploratory plots, it seemed only c3 had any noticeable effect, so you might want to concentrate on it. You would plot ty2 with something like this:

% Plotting ‘y(1)’ and ‘c3’ from ‘ty2’
[X,Y] = meshgrid(tspan,c3);
Cplot = squeeze(ty2(1,1,:,:,2));
figure(1)
surf(X, Y, Cplot)
xlabel('Time')
ylabel('C_3')
zlabel('Y_1')
title('Y_1 as a function of C_3 and Time')
grid on

You may also want to look at your system over a longer time. To do that, change the second argument in tspan. So if you wanted to increase the time from 0 to 1 (where it is now) to 0 to 5, tspan becomes:

tspan = linspace(0, 5, 25);

Have fun exploring your system’s behaviour!

[Learner]

Many many thanks for your help and support. :)

It is working in a same way as I wanted, thanks again.

[Star Strider]

My pleasure!