Solving a system of higher order differential equations with ode45 ...error

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

[Karan Gill]

If I may ask, what did you refer to when writing your code?

[Eli]

Hello Karan, I found the information for solving higher order differential equations in the mathworks 'solving a system of differential equations' function board. The main useful commands for this type of problem is odetoVectorField, matlabfunction, and ode45. Here is a link to the page https://www.mathworks.com/help/symbolic/odetovectorfield.html

[Karan Gill]

Thanks! Glad to know people read my doc :) This page is also helpful: https://www.mathworks.com/help/symbolic/examples/solve-a-second-order-differential-equation-numerically.html

Lastly, if you have trouble with inputs, try reading the input arguments info in doc. Example: Here's ode45 input info.

Answer 1

[Star Strider]

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!

[Eli]

Star Strider, I saw this answer last night immediately after you posted it. The code runs beautifully and I am very thankful for your input. The Vsubs output clarifies my confusion about which initial conditions in the vector apply to which variable. Now I have to 'guess' values for Dg and D2f until g and Df both equal close to 1. I am thinking of writing a for/while/if loop for it but I believe that might over complicate things for me and I might just be better off guessing! Again thank you so much for your input it really is amazing how helpful the mathworks community is. Ill keep you updated on the progress of the solution! Thanks!

[Star Strider]

As always, my pleasure!

Before you create a loop to tune your initial conditions, consider the approaches in Optimizing a Simulation or Ordinary Differential Equation (link). I have no experience with this, although I have done curve-fitting with a differential equation system that also treated the initial conditions as parameters to be optimized. (That is in Monod kinetics and curve fitting (link). It is similar to your problem, although not directly applicable to it. I suggest it as a way to code your problem.)

[Eli]

Just checking back in, got the problem solved and was able to get the plots I was looking for. Again thank you for your input I have seen your answers and used information you have provided many times over the years so thank you very much!

[Star Strider]

Thank you for your compliment! I do my best!

As always, my pleasure!