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?

Question

Alireza 2022-3-27

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1681734-i-have-a-system-of-ode-equations-in-the-function-file-when-i-had-to-define-zero-array-i-get-an-er

评论： Alireza 2022-3-27

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

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Torsten 2022-3-27

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1681734-i-have-a-system-of-ode-equations-in-the-function-file-when-i-had-to-define-zero-array-i-get-an-er#answer_927984

在 MATLAB Online 中打开

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.

7 个评论
显示 5更早的评论隐藏 5更早的评论

Torsten 2022-3-27

在 MATLAB Online 中打开

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