ODE solving ERROR with 5 eq
1 次查看(过去 30 天)
显示 更早的评论
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 个评论
回答(1 个)
Star Strider
2019-8-8
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 个评论
另请参阅
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!