I am trying to solve these coupled equations using ode45...I have attached both the function file and script file...

function yp=rockyruban(t,y)

a=0.5179; b=1.1003e-3; c=0.4576; d=1; hstress=3463.46; hstar=0.7846; kc=0.1137; T=503; sigma=180;

yp(1) = a*sinh((b*(sigma^d)*(1-y(2)))/((1-y(3))*(1-y(4)))) ; yp(2)= ((hstress*y(1))/(sigma^d))*(1-(y(2)/hstar)); yp(3)= (kc/3)*power((1-y(3)),4); yp(4)= c*y(1);

t0=0; tf=48; y0=[0; 0; 0; 0]; [t,y]=ode45(@rockyruban,[t0 tf],y0); plot(t,y(1)) grid xlabel('Time(hrs)') ylabel('Creepstrain(%)') title('Creep curve') hold on

Error using odearguments (line 90) ROCKYRUBAN must return a column vector.

Error in ode45 (line 113) [neq, tspan, ntspan, next, t0, tfinal, tdir, y0, f0, odeArgs, odeFcn, ...

Error in Creepruban (line 4) [t,y]=ode45(@rockyruban,[t0 tf],y0);

Answer 1

[David Goodmanson]

Hi ruban,

this addresses the immediate question (only). In the function, just before you start defining yp(1) = ..., etc, if you add the command

yp = zeros(4,1)

then the code runs. (To run this as a script I also put the function definition at the bottom, with an 'end' command as the last line of the function).

Without the preallocation, when Matlab builds up an undefined vector with x(1) = whatever, x(2) = whatever, and so forth, it assumes a row vector by default.

Unfortunately, now that it runs, ode45 can't get it done here. It issues a warning and quits. But that's a separate issue.

Answer 2

[Stephan]

Hi,

write your function in a way that it returns a column-vector like below. Then fix the error regarding your plot, which appears after fixing the first error. After this fixing the error that you built in with indexing of y in the plot-function.

Then you get this code:

t0=0;
tf=48;
y0=[0; 0; 0; 0]; 
[t,y]=ode45(@rockyruban,[t0 tf],y0);
plot(t,y(:,1))
ax = gca;
grid
xlabel('Time(hrs)')
ylabel('Creepstrain(%)')
title('Creep curve')
hold on

function yp=rockyruban(t,y)
a=0.5179;
b=1.1003e-3;
c=0.4576; 
d=1; 
hstress=3463.46;
hstar=0.7846;
kc=0.1137;
T=503;
sigma=180;
yp = [a*sinh((b*(sigma^d)*(1-y(2)))/((1-y(3))*(1-y(4))));...
    ((hstress*y(1))/(sigma^d))*(1-(y(2)/hstar));...
    (kc/3)*power((1-y(3)),4);...
    c*y(1)];
end

Problem now - The code is running, but you get a warning:

Warning: Failure at t=9.523543e+00.  Unable to meet integration
tolerances without reducing the step size below the smallest value
allowed (2.842171e-14) at time t. 

> In ode45 (line 360)
  In 'your_script_name' (line 4)

In your function code there are two warnings from matlab and i think the function has at least one content error in it:

1. You define T=503 - but you dont use T in your function. Why?
2. Your function is declared as depending from y and t - But you do not use t in the function - Why?

I could imagine that for this reason you get the warning at t=9.523543e+00 shown above, since you are not doing anything with t in your function. But this is only a guess.

However, when you change tf=48 to tf=9 you get:

.

When changing t0 = 9 and tf = 48 the warning appears at t=1.852354e+01 and your y-values explode:

.

I suggest to check the function (or maybe this is a normal bahavior for your Problem(?)) - the syntax errors that appeared are fixed as far and it is running. Regarding the warning which appears by higher tf-values i have not enough insight in the structure of your problem.

Best regards

Stephan