ODE solving ERROR with 5 eq

Question

MartinM 2019-8-8

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/475438-ode-solving-error-with-5-eq

评论： Star Strider 2019-8-12

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?

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Star Strider 2019-8-8

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/475438-ode-solving-error-with-5-eq#answer_386827

在 MATLAB Online 中打开

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')

8 个评论
显示 6更早的评论隐藏 6更早的评论

MartinM 2019-8-9

在 MATLAB Online 中打开

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 2019-8-9

在 MATLAB Online 中打开

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 2019-8-11

在 MATLAB Online 中打开

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 2019-8-11

在 MATLAB Online 中打开

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.