ODE solving ERROR with 5 eq

Hi,

I can't understant the error message from a long but simple code....

clc
clear all
close all
format long

c=3e8;
r=30e-6;
Aeff=pi.*r.*r;
E=1e-6;
lambda0=1030e-9;
om0=2.*c./lambda0;
C0=0;
eta=1E15; %Value can be changed
b2=4e-28;
b3=1.5e-45;


T0=100./om0;
Tr=20./om0;

n2=1e-23;

gamma=2.*pi.*n2./(lambda0.*Aeff);

syms Tp(z) C(z) T(z) Om(z) phi(z)

ode1 = diff(Tp) == (b2-b3.*Om).*C./Tp
ode2 = diff(C) == (b2-b3.*Om).*( (1+C.*C)./Tp.^2 ) + gamma.*E./(sqrt(2.*pi).*Tp) .*(1- Om./om0 )
ode3 = diff(T) == -b2.*Om + b3./2 .*(Om.*Om + (1+C.*C)./2.*Tp.^2 ) + 3.*gamma.*E./(2.*sqrt(2.*pi).*om0.*Tp )
ode4 = diff(Om) == gamma.*E./(sqrt(2.*pi).*Tp.^3) .*(Tr-C./om0) - eta.*E./(sqrt(2.*pi).*Tp)
ode5 = diff(phi) == 0.5.*b2.*(1./(Tp).^2 - Om.^2) + b3.*Om./3 .*(Om.*Om + 3./4 .* (C.^2 -1)./ Tp.^2 ) + 3.*gamma.*E.*(1+Om./om0) ./(4.*sqrt(2.*pi).*Tp) - 0.5.*eta.*E;

odes = [ode1; ode2; ode3; ode4; ode5]


cond1 = Tp(0) == T0;
cond2 = C(0) == 0;
cond3 = T(0) == 0;
cond4 = Om(0) == 0;
cond5 = phi(0) == 0;

conds = [cond1; cond2; cond3; cond4; cond5];
[TpSol(z), CSol(z), TSol(z), OmSol(z), phiSol(z)] = dsolve(odes,conds)

error message =

Error using sym/subsindex (line
685)
Invalid indexing or function
definition. When defining a
function, ensure that the body of
the function is a SYM object. When
indexing, the input must be
numeric, logical or ':'.

Error in ODE (line 45)
[TpSol(z), CSol(z), TSol(z),
OmSol(z), phiSol(z)] =
dsolve(odes,conds)

Do you have an idea?

Answer 1

[Star Strider]

The error was with:

[TpSol(z), CSol(z), TSol(z), OmSol(z), phiSol(z)] = dsolve(odes,conds)

since MATLAB assumed that ‘TpSol(z)’ and the rest were either function calls or that ‘z’ is an index.

However your system does not have an analytical solution. You will have to integrate it numerically.

Try this:

[VF,Sbs] = odeToVectorField(odes)
odefcn = matlabFunction(VF, 'Vars',{T,Y})

tspan = linspace(0, 100);
ics = zeros(1,4)+eps;
[t,y] = ode15s(odefcn, tspan, ics);

figure
plot(t, y)
grid
lgndc = sprintfc('%s', Sbs);
legend(lgndc, 'Location','E')

[MartinM]

Thanks for you answer.

But I am not sure to understand the sequence. I have read otehr post but can't write correctly the program.

I need to add 2 variable?

https://fr.mathworks.com/matlabcentral/answers/421503-explicit-answer-error-with-dsolve

https://fr.mathworks.com/matlabcentral/answers/426153-how-to-solve-this-ode

I add Y and K, but not working...

Not sure to understand the idea.

Many thanks again

clc
clear all
close all
format long

c=3e8;
r=30e-6;
Aeff=pi.*r.*r;
E=1e-6;
lambda0=1030e-9;
om0=2.*c./lambda0;
C0=0;
eta=1E15;
b2=-4e-28;
b2=-b2;
b3=1.5e-45;
% Tr=0.5e-15;
g=sign(b2);

T0=100./om0;
Tr=20./om0;

n2=1e-23;

gamma=2.*pi.*n2./(lambda0.*Aeff);

syms Tp(z) C(z) T(z) Om(z) phi(z) Y K

ode1 = diff(Tp) == (b2-b3.*Om).*C./Tp
ode2 = diff(C) == (b2-b3.*Om).*( (1+C.*C)./Tp.^2 ) + gamma.*E./(sqrt(2.*pi).*Tp) .*(1- Om./om0 )
ode3 = diff(T) == -b2.*Om + b3./2 .*(Om.*Om + (1+C.*C)./2.*Tp.^2 ) + 3.*gamma.*E./(2.*sqrt(2.*pi).*om0.*Tp )
ode4 = diff(Om) == gamma.*E./(sqrt(2.*pi).*Tp.^3) .*(Tr-C./om0) - eta.*E./(sqrt(2.*pi).*Tp)
ode5 = diff(phi) == 0.5.*b2.*(1./(Tp).^2 - Om.^2) + b3.*Om./3 .*(Om.*Om + 3./4 .* (C.^2 -1)./ Tp.^2 ) + 3.*gamma.*E.*(1+Om./om0) ./(4.*sqrt(2.*pi).*Tp) - 0.5.*eta.*E;

odes = [ode1; ode2; ode3; ode4; ode5]

cond1 = Tp(0) == T0;
cond2 = C(0) == 0;
cond3 = T(0) == 0;
cond4 = Om(0) == 0;
cond5 = phi(0) == 0;
conds = [cond1; cond2; cond3; cond4; cond5];

[VF,Sbs] = odeToVectorField(odes)
odefcn = matlabFunction(VF, 'Vars',{z,Y,K})
tspan = linspace(0, 100);
ics = zeros(1,4)+eps;
[Y,K] = ode15s(odefcn, tspan, ics);
figure
plot(t, y)
grid
lgndc = sprintfc('%s', Sbs);
legend(lgndc, 'Location','E')

[Star Strider]

My code worked when I posted it.

You changed the differential equation system and some of the variables in parts of the code, however you did not change other parts of my code to be compatible with those changes. I also do not understand what you are doing with the plot call. The ODE integrators output the independent variable as the first vector (usually time), and the matrix of integrated differential equations as the second argument, so the plot calls need to use the same order if you are to understand the output.

Also, while it is acceptable to pass other arguments to the ODE functions if you need to, the ODE integrators require the the first argument to the ODE function be the independent variable, and the second argument the vector of values to be integrated, so usually it would appear as ‘@(T,Y)’ (if you are integrating ’Y’ with respect to ‘T’). In your code, the ODE function is constructed as:

odefcn =

  function_handle with value:

    @(z,Y,K)[((Y(4).*1.716666666666667e-15-1.0).* ...

so ode15s is going to consider ‘z’ your independent variable, integrate ‘Y’ with respect to ‘z’, and completely ignore ‘K’ (that doesn’t appear in any of your equations anyway). The initial conditions (‘conds’, all essentiall zero) appear in the ode15s code as ‘ics’ vector.

I have no idea what you are doing with your differential equation system, however this code (with your new system of differential equations and appropriate changes) works:

c=3e8;
r=30e-6;
Aeff=pi.*r.*r;
E=1e-6;
lambda0=1030e-9;
om0=2.*c./lambda0;
C0=0;
eta=1E15;
b2=-4e-28;
b2=-b2;
b3=1.5e-45;
% Tr=0.5e-15;
g=sign(b2);
T0=100./om0;
Tr=20./om0;
n2=1e-23;
gamma=2.*pi.*n2./(lambda0.*Aeff);
syms Tp(z) C(z) T(z) Om(z) phi(z) Y  
ode1 = diff(Tp) == (b2-b3.*Om).*C./Tp
ode2 = diff(C) == (b2-b3.*Om).*( (1+C.*C)./Tp.^2 ) + gamma.*E./(sqrt(2.*pi).*Tp) .*(1- Om./om0 )
ode3 = diff(T) == -b2.*Om + b3./2 .*(Om.*Om + (1+C.*C)./2.*Tp.^2 ) + 3.*gamma.*E./(2.*sqrt(2.*pi).*om0.*Tp )
ode4 = diff(Om) == gamma.*E./(sqrt(2.*pi).*Tp.^3) .*(Tr-C./om0) - eta.*E./(sqrt(2.*pi).*Tp)
ode5 = diff(phi) == 0.5.*b2.*(1./(Tp).^2 - Om.^2) + b3.*Om./3 .*(Om.*Om + 3./4 .* (C.^2 -1)./ Tp.^2 ) + 3.*gamma.*E.*(1+Om./om0) ./(4.*sqrt(2.*pi).*Tp) - 0.5.*eta.*E;
odes = [ode1; ode2; ode3; ode4; ode5]

[VF,Sbs] = odeToVectorField(odes)
odefcn = matlabFunction(VF, 'Vars',{z,Y})
tspan = linspace(0, 10);
ics = zeros(1,5)+eps;
[Z,Y] = ode15s(odefcn, tspan, ics);
figure
plot(Z, Y)
grid
lgndc = sprintfc('%s', Sbs);
legend(lgndc, 'Location','E')

It may never produce the result you want, because the values for many of the constants are beyond the range of effective floating-point precision, in calculations involving such different magnitudes as ‘b3’ and ‘eta’.

[Star Strider]

MartinM’s Answer moved here:

Hi,

I did'nt explain the aim . This system of ODE is the linked to the propagation of a pulse (electric filed) in an optical fiber. Linked to laser-matter interrection; dispersion (same as prism) b2, interaction with gas (geeration of plasma, eta) etc.All this parameter are linked together. This come from a scientific paper, russian one, can send you if you are interrested.

Sorry, I have change stuff because the First code you send was not working for me, so I try to found other topic related to this, and try to change. . The First version you send to me was this one :, and not working

clc
clear all
close all
format long
c=3e8;
r=30e-6;
Aeff=pi.*r.*r;
E=1e-6;
lambda0=1030e-9;
om0=2.*c./lambda0;
C0=0;
eta=1E15;
b2=-4e-28;
b2=-b2;
b3=1.5e-45;
% Tr=0.5e-15;
g=sign(b2);
T0=100./om0;
Tr=20./om0;
n2=1e-23;
gamma=2.*pi.*n2./(lambda0.*Aeff);

syms Tp(z) C(z) T(z) Om(z) phi(z) 

ode1 = diff(Tp) == (b2-b3.*Om).*C./Tp
ode2 = diff(C) == (b2-b3.*Om).*( (1+C.*C)./Tp.^2 ) + gamma.*E./(sqrt(2.*pi).*Tp) .*(1- Om./om0 )
ode3 = diff(T) == -b2.*Om + b3./2 .*(Om.*Om + (1+C.*C)./2.*Tp.^2 ) + 3.*gamma.*E./(2.*sqrt(2.*pi).*om0.*Tp )
ode4 = diff(Om) == gamma.*E./(sqrt(2.*pi).*Tp.^3) .*(Tr-C./om0) - eta.*E./(sqrt(2.*pi).*Tp)
ode5 = diff(phi) == 0.5.*b2.*(1./(Tp).^2 - Om.^2) + b3.*Om./3 .*(Om.*Om + 3./4 .* (C.^2 -1)./ Tp.^2 ) + 3.*gamma.*E.*(1+Om./om0) ./(4.*sqrt(2.*pi).*Tp) - 0.5.*eta.*E;
odes = [ode1; ode2; ode3; ode4; ode5]

cond1 = Tp(0) == T0;
cond2 = C(0) == 0;
cond3 = T(0) == 0;
cond4 = Om(0) == 0;
cond5 = phi(0) == 0;
conds = [cond1; cond2; cond3; cond4; cond5];

[VF,Sbs] = odeToVectorField(odes)
odefcn = matlabFunction(VF, 'Vars',{T,Y})
tspan = linspace(0, 100);
ics = zeros(1,4)+eps;
[t,y] = ode15s(odefcn, tspan, ics);
figure
plot(t, y)
grid
lgndc = sprintfc('%s', Sbs);
legend(lgndc, 'Location','E')

And Y was not a variable, so I try to change Y, and K etc, forget my idea where I put Y and K, it was just a wrong try.

The Last code you send is working, many thanks., But Is Y(:,1) = Tp, Y(:,2)=C etc ?

The value look far from something physicaly right, is this due to my initial value b2 b3 eta etc?

I am very far to understand how it work, but I am currious. Could you explain this lines?....

[VF,Sbs] = odeToVectorField(odes)

Creats table, with equation : VF, and Name of it Sbs?

odefcn = matlabFunction(VF, 'Vars',{z,Y})

Give to Matlab the equation, and variable to be used?

tspan = linspace(0, 10);

My space of integration?

ics = zeros(1,5)+eps;

What is it for?

[Z,Y] = ode15s(odefcn, tspan, ics);

Solve everything?

MAny thanks again for you time and explanation

[Star Strider]

My pleasure.

Every Answer or Comment I posted had working code! You then changed parts of my code, so of course it then did not work. Please settle on one system of differential equations so I don’t have to chase a moving target.

This is the full version of the first code I posted, and it works:

c=3e8;
r=30e-6;
Aeff=pi.*r.*r;
E=1e-6;
lambda0=1030e-9;
om0=2.*c./lambda0;
C0=0;
eta=1E15;
b2=-4e-28;
b2=-b2;
b3=1.5e-45;
% Tr=0.5e-15;
g=sign(b2);
T0=100./om0;
Tr=20./om0;
n2=1e-23;
gamma=2.*pi.*n2./(lambda0.*Aeff);
syms Tp(z) C(z) T(z) Om(z) phi(z) T Y 
ode1 = diff(Tp) == (b2-b3.*Om).*C./Tp
ode2 = diff(C) == (b2-b3.*Om).*( (1+C.*C)./Tp.^2 ) + gamma.*E./(sqrt(2.*pi).*Tp) .*(1- Om./om0 )
ode3 = diff(T) == -b2.*Om + b3./2 .*(Om.*Om + (1+C.*C)./2.*Tp.^2 ) + 3.*gamma.*E./(2.*sqrt(2.*pi).*om0.*Tp )
ode4 = diff(Om) == gamma.*E./(sqrt(2.*pi).*Tp.^3) .*(Tr-C./om0) - eta.*E./(sqrt(2.*pi).*Tp)
ode5 = diff(phi) == 0.5.*b2.*(1./(Tp).^2 - Om.^2) + b3.*Om./3 .*(Om.*Om + 3./4 .* (C.^2 -1)./ Tp.^2 ) + 3.*gamma.*E.*(1+Om./om0) ./(4.*sqrt(2.*pi).*Tp) - 0.5.*eta.*E;
odes = [ode1; ode2; ode3; ode4; ode5]
% cond1 = Tp(0) == T0;
% cond2 = C(0) == 0;
% cond3 = T(0) == 0;
% cond4 = Om(0) == 0;
% cond5 = phi(0) == 0;
% conds = [cond1; cond2; cond3; cond4; cond5];
[VF,Sbs] = odeToVectorField(odes)
odefcn = matlabFunction(VF, 'Vars',{T,Y})
tspan = linspace(0, 100);
ics = zeros(1,4)+eps;
[t,y] = ode15s(odefcn, tspan, ics);
figure
plot(t, y)
grid
lgndc = sprintfc('%s', Sbs);
legend(lgndc, 'Location','E')

With respect to the functions, please refer to the documentation for a full description of what they do and how they work.

First:

[VF,Sbs] = odeToVectorField(odes)

This takes your system of ODEs and converts it to a system of first-degree ODEs so that the matlabFunction function can convert it to an anonymous function (here) that the numerical ODE solvers can work with. Note that odeToVectorField uses ‘Y’ as its dependent variable, so use the second output (here ‘Sbs’) to understand the substitutions it made, and what ‘Y[1]’, ‘Y[2]’ and the others correspond to in the original differential equation system.

See odeToVectorField for details.

Second:

odefcn = matlabFunction(VF, 'Vars',{z,Y})

This creates an anonymous function from its arguments (‘VF’ here) that MATLAB ca use with the numeric ODE solvers. Neither odeToVectorField nor matlabFunction assumes that an independent variable exists, so you need to provide one in the 'Vars' argument so the ODE integrators can use the function correctly.

See: matlabFunction and Anonymous Functions for details.

Third:

tspan = linspace(0, 10);

This creates a vector of independent variable elements (100 by default) that the ODE solvers use to define the limits of integration.

See linspace for details.

Fourth:

ics = zeros(1,5)+eps;

This defines the initial conditions vector (essentially ‘conds’ but with numeric, not symbolic, elements).

Fifth:

[Z,Y] = ode15s(odefcn, tspan, ics);

This integrates the system of first-degree ODEs created by odeToVectorField and matlabFunction. The independent variable the function creates is ‘Z’ (here), and the integrated equations are ‘Y’ (here). The ode15s function is particularly useful for ‘stiff’ systems, such as yours, where the constants vary by several orders-of-magnitude. Note that you can call the outputs anything you want. I chose ‘Z’ and ‘Y’.)

See: ode15s for details.

[Star Strider]

MartinM’s Answer moved here:

Please use ‘Comment on this Answer’ — not ‘Answer this question’ — for comments.

—————

That's clear

But the code can't be run in my computer. It stoped at

odefcn = matlabFunction(VF, 'Vars',{T,Y})

with error

Error using
sym/matlabFunction>checkVars
(line 159)
Variable names must be valid
MATLAB variable names.

Error in sym/matlabFunction (line 104)
vars = checkVars(funvars,opts);

If I replace T by z it run...?

OTHER question, Why my initials conditionare not usefull. For me they should be the key to resolve the system...

Regard

[Star Strider]

The code I posted runs without error in R2019a.

What variable names are you using?

Did you change something again? You did not post the code you used that threw that error, so I cannot help you with that specific problem.

‘OTHER question, Why my initials condition are not usefull.’

Your ‘conds’ variable is symbolic and contains symbolic values. The numeric ODE solvers cannot use symbolic variables or values, so I created a numeric vector for the initial conditions that ode15s can use.

[Torsten]

@Star Strider

Defining T(z) and T as syms might cause a problem ?

[Star Strider]

@Torsten — Thank you! It definitely could. I didn’t see that (early here).