when solving Gilmore equation using ode45 , the answer comes as zero and it plots a straight line in the time interval

clc

clear all

close all

%

q_0 = 998;

C_0 = 1484; %changes did 8000

r = 6.5; %changes did 6.5

mu = 0.001;

sigma = 0.072;

p_0 = 101;

R_0 = 61;%0.000061;

m = 0.11;

k = 0;

b = 1/(q_0 * (C_0)^2);

B = ((b*r - p_0))^(-1);

syms R(t)

eqn = (1- (diff(R(t),t)/(sqrt(C_0^2 + (r-1)*(((b*r)^(-1/r))/q_0)*(r/r-1)*(((p_0 + 2*sigma/R_0)* ((R_0^3 - m* R_0^3)/(R^3 - m* R_0^3))^k - (4*mu*diff(R,1)/R) - 2*sigma/R - B)^(1- 1/r) - (p_0 + B)^(1- 1/r))))))*R* diff(R(t),t,t) + (3/2)*(1- (diff(R(t),t)/(3*(sqrt(C_0^2 + (r-1)*(((b*r)^(-1/r))/q_0)*(r/r-1)*(((p_0 + 2*sigma/R_0)* ((R_0^3 - m* R_0^3)/(R^3 - m* R_0^3))^k - (4*mu*diff(R,1)/R) - 2*sigma/R - B)^(1- 1/r) - (p_0 + B)^(1- 1/r)))))))*((diff(R(t),t))^2) == (1+ (diff(R(t),t)/(sqrt(C_0^2 + (r-1)*(((b*r)^(-1/r))/q_0)*(r/r-1)*(((p_0 + 2*sigma/R_0)* ((R_0^3 - m* R_0^3)/(R^3 - m* R_0^3))^k - (4*mu*diff(R,1)/R) - 2*sigma/R - B)^(1- 1/r) - (p_0 + B)^(1- 1/r))))))*((((b*r)^(-1/r))/q_0)*(r/r-1)*(((p_0 + 2*sigma/R_0)* ((R_0^3 - m* R_0^3)/(R^3 - m* R_0^3))^k - (4*mu*diff(R,1)/R) - 2*sigma/R - B)^(1- 1/r) - (p_0 + B)^(1- 1/r))) + (1- (diff(R(t),t)/(sqrt(C_0^2 + (r-1)*(((b*r)^(-1/r))/q_0)*(r/r-1)*(((p_0 + 2*sigma/R_0)* ((R_0^3 - m* R_0^3)/(R^3 - m* R_0^3))^k - (4*mu*diff(R,1)/R) - 2*sigma/R - B)^(1- 1/r) - (p_0 + B)^(1- 1/r))))))*(R/(sqrt(C_0^2 + (r-1)*(((b*r)^(-1/r))/q_0)*(r/r-1)*(((p_0 + 2*sigma/R_0)* ((R_0^3 - m* R_0^3)/(R^3 - m* R_0^3))^k - (4*mu*diff(R,1)/R) - 2*sigma/R - B)^(1- 1/r) - (p_0 + B)^(1- 1/r)))))* (((b*r)^(-1/r))/q_0 * ((((p_0 + 2*sigma/R_0)* ((R_0^3 - m* R_0^3)/(R^3 - m*R_0^3))^k - (4*mu*diff(R,1)/R) - 2*sigma/R - B)^(-1/r))*((p_0 + 2*sigma/R_0)*k*(((R_0^3 - m* R_0^3)/(R^3 - m* R_0^3))^(k-1))*((-(R_0^3 - m*(R_0^3))*3*(R^2)*diff(R(t),t))/(R^3 - m*R_0^3)^2) - ((4*mu*R*diff(R(t),t,t) - 4*mu*(diff(R(t),t))^2)/R^2) + (2*sigma*diff(R(t),t))/R^2)));

[eqs,vars] = reduceDifferentialOrder(eqn,R(t))

initConditions = [364 0]

[M,F] = massMatrixForm(eqs,vars)

f = M\F

odefun = odeFunction(f,vars);

tspan = 000000:1:7000000;

[t,yexact] = ode15s(odefun, tspan, initConditions);

plot(t,yexact(:,1))

Answer 1

[darova]

Here is what i tried. Do you know the result should looks like?

[kirankumari sheshma]

Thank you very much. That is a much better result than what i have tried to get. The result of the graph should look like in the image.png that i have uploaded.

[kirankumari sheshma]

if you go to page no.147 of the uploaded pdf fille, the problem is there.

Thank you in advance.

[darova]

Here is a better result

What about complex results? Take look also on units

>> (0.1-9.9)^((r-1)/r)
ans =
  -6.1080 + 3.2057i

[kirankumari sheshma]

If you are changing KPa to Pa , i think so it will be multiplied by 1000, instead of dividing by thousand. ?

The complex number will not be an issue, i guess.

The graph is supposed to rebound after 65 microseconds. Do you have any idea as to why the grpah is not rebounding ?

[darova]

• If you are changing KPa to Pa , i think so it will be multiplied by 1000, instead of dividing by thousand. ?

Absolutely corrrect

P0 = 101e3;     % [Pa]
plot(t,abs( u(:,1)*1e6 ))

Looks like bounce, but not at

[kirankumari sheshma]

I think so, that the simulation done in the paper along with the data do not match the given idea. I doubt on the results being published in the paper.

Thank you very much for your help.

[KALYAN ACHARJYA]

Can you refer the paper link?

[kirankumari sheshma]

the uploaded PDF is the given paper.

Answer 2

[KALYAN ACHARJYA]

"it plots a straight line in the time interval"

Because yexact(:,1) values return as same value for all 1st column elemnets

>>yexact(:,1) 
   364
   364
   364
   364
   364
   364
   ....
   ...

As t values are changing across x axes, wheare as yexact(:,1) data elements are same

[kirankumari sheshma]

That is my question. I would like to know why the yexact values are same ? when i use ode45 function, it still comes as the same.

i have given a picture image of the gilmore equation problem that i want to solve .