MATLAB: Solving a system of higher order differential equations with ode45 …error

differential equationsode45

Hello. I am an engineering student working on a rather simple problem but have exhausted my searches online and am not sure what is wrong with my setup. I have two equations, one third order and one second order differential equation. I also have 5 initial conditions where two of the initial conditions are guesses (namely f''(0) and g'(0)) – (typical guess between .5-1) when these guesses are correct the values of f'(large n) and g(large n) will converge to 1. The other three IC's are f(0)=0, g(0)=2.5, f'(0)=0. I wrote a code to reduce the order of the differential equations to a system of first order differential equations using the odeToVectorField command, then using matlab function and ode45 to 'solve' the system simultaneously… I likely have more than one mistake in here but am at a loss…I am getting NAN, complex doubles, and other nonsense outputs for the Sol structure…Any help would be appreciated. Also please let me know if I need to provide more information. Thank you.
    % code
clear all
close all
syms g(n) f(n)
%
C=g^-(1/3);
Pr=.7;
gamma=1.4;
Me=3;
gw=1+sqrt(Pr)*((gamma-1)/2)*(Me^2) %Taw/T_e
EQ1=diff(C*diff(f,2))==-f*diff(f,2);
EQ2=diff((C/Pr)*diff(g))+f*diff(g)+(C*(Me^2*(gamma*(gamma-1))))*diff(f,2)^2==0;
EQS=[EQ1 EQ2];
V = odeToVectorField(EQS)
n=[0:1:10^2];
M = matlabFunction(V,'vars', {'n','Y'});
ICs=[f(0)==0,diff(f(0))==0,diff(f(0),2)==.75,g(0)==gw,diff(g(0))==.75] %not sure how to apply the 5 IC's to correspond to V - when i run this i get error div/0 - but oddly enough if I change around the order of the IC's the error is replaced with inputs must be type double, but then I get all NAN results I must be doing something really wrong. %
Sol = ode45(M,[n(1) n(end)],ICs)
  end

Best Answer

Try this:
syms g(n) f(n) n Y 
%
C=g^-(1/3);
Pr=.7;
gamma=1.4;
Me=3;
gw=1+sqrt(Pr)*((gamma-1)/2)*(Me^2) %Taw/T_e
EQ1=diff(C*diff(f,2))==-f*diff(f,2);
EQ2=diff((C/Pr)*diff(g))+f*diff(g)+(C*(Me^2*(gamma*(gamma-1))))*diff(f,2)^2==0;
EQS=[EQ1 EQ2];
[V, Vsubs] = odeToVectorField(EQS)
n=[0:1:1E+2];
M = matlabFunction(V,'vars', {'n','Y'});
% ICs=[f(0)==0,diff(f(0))==0,diff(f(0),2)==.75,g(0)==gw,diff(g(0))==.75] %not sure how to apply the 5 IC's to correspond to V - when i run this i get error div/0 - but oddly enough if I change around the order of the IC's the error is replaced with inputs must be type double, but then I get all NAN results I must be doing something really wrong. %
ICs = [gw; 0.75; 0; 0.75; 0.75];
Sol = ode45(M,[n(1) n(end)],ICs)
figure(1)
semilogy(Sol.x, abs(Sol.y))
This runs, although you will probably need to be certain all the variables have the appropriate initial conditions.
I tweaked your code just a bit, adding ‘Vsubs’ to the odeToVectorField output to guarantee that I understand the ‘Y’ variables that correspond to the functions and derivatives. I then did the best I could to match them with the ‘ICs’ vector.
I added the plot. You will need to change it to provide the correct (and easily understandable) variables.
I didn’t see this earlier today, so I thank you for bringing it to my attention. (I generally don’t reply to private messages.) Your Question is straightforward, well-stated, and (I hope) easily solved!
Related Question