I have a system of ODE equations. In the function file when I had to define Zero Array, I get an error for matrix array. How can I fix it?

clc;
clear;
close all;
global k1 k2 k3 k4 Cp Re k We De MuR delt n bta yii www

Cp=0.5;
Re=.5;
We=.2;
De=1;
k1=10;  
k2=1;     
k3=.5;
k4=1;

delt=2;
bta=0.37;
k=1.4;
n=180;

a=zeros(n,1); b=zeros(n,1); a(1)=1; b(1)=1;
for k=2:n, a(k)=2*(k-1)+1; b(k)=(k-1)^2; end
JacM=diag(a)+diag(sqrt(b(2:n)),1)+diag(sqrt(b(2:n)),-1);
[w,x]=eig(JacM); x=diag(x); w=w(1,:)'.^2;
yii=x; www=w;

MuRi=[0.1];
np=numel(MuRi);
figure(1);
COLORS=hsv(np);
for i=1:np
    MuR=MuRi(i);
    options = odeset('RelTol',1e-9,'AbsTol',1e-9,'MaxStep',1e-9);
    iniMat=[1,0,ones(1,n)];
    [time,Y] = ode15s(@FthixoPrf,[0 80],iniMat);
    
    subplot(2,1,1);
    plot(time,Y(:,1),'Color',COLORS(i,:),'DisplayName',['MuR = ' num2str(MuR)]);
    hold on;
    subplot(2,1,2);
    plot(time,Y(:,3),'Color',COLORS(i,:));
    hold on;
end
grid on;

unction file:

function dy=FthixoPrf(t,y)
    global k1 k2 k3 k4 Cp Re k We De MuR delt n bta yii www
    dy = zeros(n,n);
    R=y(1);U=y(2);
    Str=y(3:n+2);
    gamRR=y(n+3:2*n+2);
    gamTT=y(2*n+3:3*n+2)
    dy(1)=U;
    sum=0;
    for i=1:n
         s1=www(i)*exp(yii(i))*(MuR+Str(i))/(yii(i)+R^3)^2;
         s2=www(i)*exp(yii(i))*Str(i)*(gamRR(i)-gamTT(i))/(yii(i)+R^3);
         sum=(-4*R*U*s1/Re)+(2*s2/(3*Re*De*R));
    end
    dy(2)=(-1.5*U^2+Cp*((1+We)*(1/R)^(3*k)-We/R-(1+delt*sin(t))))/R+sum;

    for i=1:n
         epsi=abs(2*sqrt(3)*R^2*U/(R^3+yii(i)));
         dy(2+i)=(1/(t+0.0001)^bta)*(-k1*Str(i)*epsi+k2*(epsi)^0.5*(1-Str(i))+k3*(1-Str(i))); 
         inv2sqrt=(k2*(epsi)^0.5*+k3)/(k1*epsi+k2*(epsi)^0.5+k3);
         dy(n+2+i)=(k4^bta/(t+0.0001)^bta)*(MuR+inv2sqrt)*gamRR(i)*(4*R^2*U/(R^3+yii(i)));
         dy(2*n+2+i)=-(k4^bta/(t+0.0001)^bta)*(MuR+inv2sqrt)*gamTT(i)*(2*R^2*U/(R^3+yii(i)));
     end
end

Answer 1

[Torsten]

dy = zeros(3*n+2,1)

But then you also have to specify 3*n+2 initial values for the variables you solve for.

You only specify n+2.

[Alireza]

After changin dy, I stil have that error for gamRR and gamTT:

Index exceeds matrix dimensions.

[Torsten]

Did you specify a vector of length 3*n+2 as iniMat ?

[Alireza]

No I didnt.

How must exactly I must do it to be true?

[Torsten]

You must have as many solution functions as you define equations for them.

I don't know which equations you try to solve and what are senseful initial conditions.

I only know that if you have 3*n+2 derivatives, you usually have also 3*n+2 functions to be determined.

[Alireza]

yes for that I wrote functions in function file: dy(2+i),dy(n+2+i), dy(2*n+2+i)

But it still has error: Index exceeds matrix dimensions.

and I dont know how to solve it.

[Torsten]

I don't know your initial conditions. I just changed them to

iniMat=[1;0;ones(3*n,1)];

and in the function

dy = zeros(3*n+2,1);

and the code is syntactically correct.

global k1 k2 k3 k4 Cp Re k We De MuR delt n bta yii www
Cp=0.5;
Re=.5;
We=.2;
De=1;
k1=10;  
k2=1;     
k3=.5;
k4=1;
delt=2;
bta=0.37;
k=1.4;
n=180;
a=zeros(n,1); b=zeros(n,1); a(1)=1; b(1)=1;
for k=2:n, a(k)=2*(k-1)+1; b(k)=(k-1)^2; end
JacM=diag(a)+diag(sqrt(b(2:n)),1)+diag(sqrt(b(2:n)),-1);
[w,x]=eig(JacM); x=diag(x); w=w(1,:)'.^2;
yii=x; www=w;
MuRi=[0.1];
np=numel(MuRi);
figure(1);
COLORS=hsv(np);
for i=1:np
    MuR=MuRi(i);
    options = odeset('RelTol',1e-9,'AbsTol',1e-9,'MaxStep',1e-9);
    iniMat=[1;0;ones(3*n,1)];
    [time,Y] = ode15s(@FthixoPrf,[0 80],iniMat);
    
    subplot(2,1,1);
    plot(time,Y(:,1),'Color',COLORS(i,:),'DisplayName',['MuR = ' num2str(MuR)]);
    hold on;
    subplot(2,1,2);
    plot(time,Y(:,3),'Color',COLORS(i,:));
    hold on;
end
grid on

function dy=FthixoPrf(t,y)

global k1 k2 k3 k4 Cp Re k We De MuR delt n bta yii www
t
dy = zeros(3*n+2,1);
R=y(1);U=y(2);
Str=y(3:n+2);
gamRR=y(n+3:2*n+2);
gamTT=y(2*n+3:3*n+2);
dy(1)=U;
sum=0;
for i=1:n
    s1=www(i)*exp(yii(i))*(MuR+Str(i))/(yii(i)+R^3)^2;
    s2=www(i)*exp(yii(i))*Str(i)*(gamRR(i)-gamTT(i))/(yii(i)+R^3);
    sum=(-4*R*U*s1/Re)+(2*s2/(3*Re*De*R));
end
dy(2)=(-1.5*U^2+Cp*((1+We)*(1/R)^(3*k)-We/R-(1+delt*sin(t))))/R+sum;
for i=1:n
    epsi=abs(2*sqrt(3)*R^2*U/(R^3+yii(i)));
    dy(2+i)=(1/(t+0.0001)^bta)*(-k1*Str(i)*epsi+k2*(epsi)^0.5*(1-Str(i))+k3*(1-Str(i))); 
    inv2sqrt=(k2*(epsi)^0.5*+k3)/(k1*epsi+k2*(epsi)^0.5+k3);
    dy(n+2+i)=(k4^bta/(t+0.0001)^bta)*(MuR+inv2sqrt)*gamRR(i)*(4*R^2*U/(R^3+yii(i)));
    dy(2*n+2+i)=-(k4^bta/(t+0.0001)^bta)*(MuR+inv2sqrt)*gamTT(i)*(2*R^2*U/(R^3+yii(i)));
end
end

[Alireza]

Thank you so much.