ODE solving ERROR with 5 eq
2 vues (au cours des 30 derniers jours)
Afficher commentaires plus anciens
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 commentaires
Réponses (1)
Star Strider
le 8 Août 2019
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 commentaires
Star Strider
le 12 Août 2019
@Torsten — Thank you! It definitely could. I didn’t see that (early here).
Voir également
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!