Index exeeded error for ode45 function

1 Ansicht (letzte 30 Tage)
AKSHAY PUJARI
AKSHAY PUJARI am 1 Dez. 2020
Beantwortet: Stephan am 1 Dez. 2020
hi all,
I have generated a ODE matrix (M)as shown in the code.
ri=0.53;
syms w(r)
Nrad(r)=-(( (r^-2) - 1)./((ri^-2)+1));
Ntheta(r)=(( (r^-2) + 1)./((ri^-2)-1));
[V] = odeToVectorField(diff(w, 4) + (2/r)*diff(w,3) - ((1/r^2) + 6*Nrad(r) )*diff(w,2) + ((((1/r^2)-6*Ntheta(r) ))/r)*diff(w,1) == 0 )
M = matlabFunction(V,'vars', {'r','Y'})
sol = ode45(M,[0.53 0.9],[0 0.53]);
fplot(@(x)deval(sol,x,1), [0.53, 0.9])
The ode45 is generating a error as follows:
what is causing this?
Index exceeds the number of array elements (2).
Error in
symengine>@(r,Y)[Y(2);Y(3);Y(4);-(1.0./r.^2.*3.157935826372082e-1-1.315793582637208).*Y(3)-(Y(4).*2.0)./r+((1.0./r.^2.*1.343763037129746+2.343763037129746).*Y(2))./r]
Error in odearguments (line 90)
f0 = feval(ode,t0,y0,args{:}); % ODE15I sets args{1} to yp0.
Error in ode45 (line 115)
odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);
Error in ODE_AnnularDisc_Nagarjuna (line 17)
sol = ode45(M,[0.53 0.9],[0 0.53]);

Akzeptierte Antwort

Stephan
Stephan am 1 Dez. 2020
This is an ode of order 4 - so odeToVectorField will reformulate it to a system of 4 ode with order 1:
>> pretty(V)
/ Y \
| 2 |
| |
| Y |
| 3 |
| |
| Y |
| 4 |
| |
| / 3221 5618 \ |
| | ------- + ---- | Y |
| | 2 2397 | 2 2 Y |
| \ 2397 r / 4 / 4045 16854 \ |
| --------------------- - ---- - | -------- - ----- | Y |
| r r | 2 12809 | 3 |
\ \ 12809 r / /
Thus you have to tell ode45 a number of 4 initial points - one for every equation. Using some fantasy values would give:
ri=0.53;
syms w(r)
Nrad(r)=-(( (r^-2) - 1)./((ri^-2)+1));
Ntheta(r)=(( (r^-2) + 1)./((ri^-2)-1));
[V] = odeToVectorField(diff(w, 4) + (2/r)*diff(w,3) - ((1/r^2) + 6*Nrad(r) )*diff(w,2) + ((((1/r^2)-6*Ntheta(r) ))/r)*diff(w,1) == 0 )
M = matlabFunction(V,'vars', {'r','Y'})
sol = ode45(M,[1 10],[0 0 0.53 0.9]);
fplot(@(x)deval(sol,x,1), [1, 10])

Weitere Antworten (0)

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by