MATLAB: How to plug the symbol vector into function.m file for ode solving

MATLABodeode45

I want to solve the problem where a stage is that I need to use ode23 or ode45 to solve the time-varying variables. For instance, the varibales is nX1 space vector. So I need to build the nX1 differential equations like that

function f = ODE(t,y,Q,R,r)

f =[2*y(1) – 10*exp(-2*t) – 22*y(3) – (y(1)*(y(1)/2 + 25*y(2) – (51*y(3))/2))/(2*(t – 15)) – (25*y(2)*(y(1)/2 + 25*y(2) – (51*y(3))/2))/(t – 15) + (51*y(3)*(y(1)/2 + 25*y(2) – (51*y(3))/2))/(2*(t – 15))

51*y(2) – 60*y(3) – 11*y(6) – (y(1)*(y(2)/2 + 25*y(5) – (51*y(6))/2))/(2*(t – 15)) – (25*y(2)*(y(2)/2 + 25*y(5) – (51*y(6))/2))/(t – 15) + (51*y(3)*(y(2)/2 + 25*y(5) – (51*y(6))/2))/(2*(t – 15))

11*y(3) – 10*y(4) – 11*y(8) – (y(1)*(y(3)/2 + 25*y(6) – (51*y(8))/2))/(2*(t – 15)) – (25*y(2)*(y(3)/2 + 25*y(6) – (51*y(8))/2))/(t – 15) + (51*y(3)*(y(3)/2 + 25*y(6) – (51*y(8))/2))/(2*(t – 15))…]

Traditionally, I can solve it by using [t,y]=ode45(@ODE,tspan,InitialCondition) to find the curve of all of the variables.

My point is: Can I let the differnential equations be automatically solved as soon as I get the symbolic expression of matrix f. I have think about it for serval days already.

Thanks.

Best Answer

Are you manually pasting your symboling expression in the ODE() function? You can use odetovectorfield(): https://www.mathworks.com/help/symbolic/odetovectorfield.html to solve the equations with ode45(). See the example "Solve Higher-Order Differential Equation Numerically" on the documentation page of odetovectorfield().

Related Solutions

MATLAB: When I run the program inaccurate results appear What is the problem

The best approach to this is likely not the vpasolve function, since it takes forever.

I have no idea of these are the ‘accurate’ results you want (I have no idea what that means), however you will get them much more quickly:

syms t y1 y2 y3 y4  yn11 yn1 yr yr11 yr111   yr v y(v) x
format longe
h=0.1;
r=1.7686999343761197980595002942562;   %1694/911=1.859495060373216e+00
x=0;
y0=0;
y01=0.25;
y011=0;
fn=-y01+y0.*y01.*y011;
fnr=-yr11+yr.*yr11.*yr111;
fn1=-yn1+y1.*yn1.*yn11;
fn2=-y3+y2.*y3.*y4;
eq1 =-y2-(r - 2)/r*y0+ 2/(r*(r - 1))*yr+ (2*r - 4)/(r - 1)*y1 -(h.^3.*(r^4 - 8*r^3 + 22*r^2 - 23*r + 6))/(120*r)*fn+ (h.^3.*(- r^2 + 2*r + 3))/(60*r)*fnr -(h.^3.*(- r^3 + 4*r^2 + 9*r - 26))/60*fn1+ (h.^3.*(- r^3 + 2*r + 1))/120*fn2;
eq1 = simplifyFraction(eq1);
eq1d = vpa(eq1, 5)
eq2= -h.*y3 -(r - 3)/r*y0+ 3/(r*(r - 1))*yr+ (r - 4)/(r - 1)*y1 -(h.^3.*(3*r^4 - 24*r^3 + 66*r^2 - 67*r + 12))/(240*r)*fn -(h.^3.*(r^4 - 5*r^3 + 5*r^2 + 5*r - 4))/(40*r*(r^2 - 3*r + 2))*fnr -(h.^3.*(- 3*r^4 + 15*r^3 + 15*r^2 - 137*r + 116))/(120*(r - 1))*fn1+ (h.^3.*(- r^4 + 2*r^3 + 2*r^2 + 3*r - 12))/(80*(r - 2))*fn2;
eq2 = simplifyFraction(eq2);
eq2d = vpa(eq2, 5)
eq3=-h.^2.*y4+2/r*y0+ 2/(r*(r - 1))*yr -2/(r - 1)*y1+ (h.^3.*(- r^4 + 8*r^3 - 22*r^2 + 18*r + 5))/(120*r)*fn -(h.^3.*(r^4 - 5*r^3 + 5*r^2 + 5*r + 5))/(60*r*(r^2 - 3*r + 2))*fnr -(h.^3.*(- r^4 + 5*r^3 + 5*r^2 - 75*r + 77))/(60*(r - 1))*fn1+ (h.^3.*(- r^4 + 2*r^3 + 2*r^2 + 42*r - 81))/(120*(r - 2))*fn2;
eq3 = simplifyFraction(eq3);
eq3d = vpa(eq3, 5)
eq4=-h.*yn1 -(r - 1)/r*y0+ 1/(r*(r - 1))*yr+ (r - 2)/(r - 1)*y1 -(h.^3.*(r^4 - 8*r^3 + 22*r^2 - 21*r + 6))/(240*r)*fn+ (h.^3.*(- r^2 + 2*r + 3))/(120*r)*fnr -(h.^3.*(- r^3 + 4*r^2 + 9*r - 8))/120*fn1 -(h.^3.*(r^3 - 2*r + 1))/240*fn2;
eq4 = simplifyFraction(eq4);
eq4d = vpa(eq4, 5)
eq5=-h.^2.*yn11+ 2/r*y0+ 2/(r*(r - 1))*yr -2/(r - 1)*y1 -(h.^3.*(r^4 - 8*r^3 + 22*r^2 - 28*r + 10))/(120*r)*fn -(h.^3.*(r^4 - 5*r^3 + 5*r^2 + 5*r - 10))/(60*r*(r^2 - 3*r + 2))*fnr -(h.^3.*(- r^4 + 5*r^3 + 5*r^2 - 35*r + 22))/(60*(r - 1))*fn1+ (h.^3.*(- r^4 + 2*r^3 + 2*r^2 - 8*r + 4))/(120*(r - 2))*fn2;
eq5 = simplifyFraction(eq5);
eq5d = vpa(eq5, 5)
eq6=-h.*yr11+ (r - 1)/r*y0+ (2*r - 1)/(r*(r - 1))*yr -r/(r - 1)*y1+ (h.^3.*(2*r^4 - 13*r^3 + 28*r^2 - 22*r + 5))/240*fn+ (h.^3.*(4*r^3 - 15*r^2 + 10*r + 5))/(120*(r - 2))*fnr+ (h.^3.*r*(- 2*r^3 + 7*r^2 + 2*r - 3))/120*fn1+ (h.^3.*r*(2*r^4 - 5*r^3 + 2*r^2 + 2*r - 1))/(240*(r - 2))*fn2;
eq6 = simplifyFraction(eq6);
eq6d = vpa(eq6, 5)
eq7=-h.^2.*yr111+ 2/r*y0+ 2/(r*(r - 1))*yr -2/(r - 1)*y1+ (h.^3.*(4*r^4 - 22*r^3 + 38*r^2 - 22*r + 5))/(120*r)*fn -(h.^3.*(- 14*r^4 + 55*r^3 - 55*r^2 + 5*r + 5))/(60*r*(r^2 - 3*r + 2))*fnr -(h.^3.*(4*r^4 - 15*r^3 + 5*r^2 + 5*r - 3))/(60*(r - 1))*fn1+ (h.^3.*(4*r^4 - 8*r^3 + 2*r^2 + 2*r - 1))/(120*(r - 2))*fn2;
eq7 = simplifyFraction(eq7);
eq7d = vpa(eq7, 5)
eq8=-h.*y01 -(r + 1)/r*y0 -1/(r*(r - 1))*yr+ r/(r - 1)*y1+ (h.^3.*(r^3 - 8*r^2 + 22*r - 5))/240*fn+ (h.^3.*(r^3 - 5*r^2 + 5*r + 5))/(120*(r^2 - 3*r + 2))*fnr+ (h.^3.*r*(- r^3 + 5*r^2 + 5*r - 3))/(120*(r - 1))*fn1 -(h.^3.*r*(- r^3 + 2*r^2 + 2*r - 1))/(240*(r - 2))*fn2;
eq8 = simplifyFraction(eq8);
eq8d = vpa(eq8, 5)
eq9 =-h.^2.*y011+ 2/r*y0+ 2/(r*(r - 1))*yr -2/(r - 1)*y1 -(h.^3.*(r^4 - 8*r^3 + 22*r^2 + 22*r - 5))/(120*r)*fn -(h.^3.*(r^4 - 5*r^3 + 5*r^2 + 5*r + 5))/(60*r*(r^2 - 3*r + 2))*fnr -(h.^3.*(- r^4 + 5*r^3 + 5*r^2 + 5*r - 3))/(60*(r - 1))*fn1+ (h.^3.*(- r^4 + 2*r^3 + 2*r^2 + 2*r - 1))/(120*(r - 2))*fn2;
eq9 = simplifyFraction(eq9);
eq9d = vpa(eq9, 5)
eqsfcn = matlabFunction([eq1,eq2,eq3,eq4,eq5,eq6,eq7,eq8,eq9], 'Vars',{[y1,y2,y3,y4,yn1,yn11,yr,yr11,yr111]})
eqsol = fsolve(eqsfcn, rand(1,9))
eqsolc = num2cell(eqsol);
[y1,y2,y3,y4,yn1,yn11,yr,yr11,yr111] = eqsolc{:}

I used the simplifyFraction function to do just that with your equations, then the fsolve function to do the actual calculations. These changes make it significantly more efficient.

The results I got are:

y1 =
     2.495835158216658e-02
y2 =
     4.966725033036157e-02
y3 =
     2.450145921764798e-01
y4 =
    -4.970810349993776e-02
yn1 =
     2.487509121241510e-01
yn11 =
    -2.496352514349295e-02
yr =
     4.398727143758653e-02
yr11 =
     2.460985502335340e-01
yr111 =
    -4.401564490668385e-02

You must be certain your equations are written as you want them to be in order to get ‘accurate’ values. We have no control over that.

Experiment to get the results you want.

MATLAB: Trying to code a tricky double summation

Something like this should work

s = 0;
for i=0:min(x,y)
    for k=0:x-1
        s = s + (-1)^k * D(2*x-2*k-2*i, 2*y+2-2*i) * T(2*x-2*k-2*i, 2*y+2-2*i)^(z/2)
    end
end

Best Answer

Related Solutions

MATLAB: When I run the program inaccurate results appear What is the problem

MATLAB: Trying to code a tricky double summation

Related Question