How to solve systems of ode in matlab?

Question

RITIKA Jaiswal 2023-4-10

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1944714-how-to-solve-systems-of-ode-in-matlab

编辑： Walter Roberson 2023-4-11

Adobe_Scan_Apr_10,_2023[1].pdf

I am trying to solve following sets of odes for different cases of control volume .I have written the code here and after running the code I am getting error .Please help me to fix the error.I have considered the value of of all K as constant in the code for now.

Lx = 0.1;
Ly = 0.1;
dx = 0.01;
dy = 0.01;
nx = Lx/dx;
ny = Ly/dy;
Tin = 500;
x = dx/2:dx:Lx+(dx/2);
y = dy/2:dy:Ly+(dy/2);
Nx = nx+1;
Ny = ny+1;
M=zeros(Nx,Ny);
Tini = 600;
W=Tini*(1+M);
Told = W ;
Counter = 0;
for i = 1:Nx
    for j = 1:Ny
        if (i==1) && (j==Ny)
            dvdt=@(t,W)[3*Tin-2*W(i,j)+4*W(i,j-1)];
            
        end
        if (i==1) && (j==1)
            dvdt=@(t,W)[3*Tin + 4*W(i,j)-3*W(i,j+1)];
        end
        if (i==Nx) && (j==1)
            dvdt=@(t,W)[3*T(i-1,j) + 4*W(i,j)-3*W(i,j+1)];
        end
        if (i==Nx) && (j==Ny)
            dvdt=@(t,W)[3*T(i-1,j-1) + 4*W(i,j)-3*W(i,j+1)];
        end
        
        if (j==Ny)
            dvdt=@(t,W)[2*W(i-1,j)+4*W(i,j-1)-3*W(i,j)];
        end
        if (j==1)
            dvdt=@(t,W)[4*W(i+1,j+1)-3*W(i,j)];
        else
            dvdt=@(t,W)[3*T(i-1,j) + 4*W(i+1,j-1)-3*W(i-1,j+1)];
        end
        [t,W]=ode45(dvdt,[0 0.4],Told)
        
    end
end
plot(t,W)

below are the Following errors which I am getting after running my code:

Attempted to access W(2,2); index out of bounds because size(W)=[121,1].

Error in @(t,W)[4*W(i+1,j+1)-3*W(i,j)]

Error in odearguments (line 87)

f0 = feval(ode,t0,y0,args{:}); % ODE15I sets args{1} to yp0.

Error in ode45 (line 113)

[neq, tspan, ntspan, next, t0, tfinal, tdir, y0, f0, odeArgs, odeFcn, ...

Error in FVM (line 41)

[t,W]=ode45(dvdt,[0 0.4],Told)

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

请先登录，再进行评论。

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

Answer 1

Torsten 2023-4-10

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1944714-how-to-solve-systems-of-ode-in-matlab#answer_1213074

在 MATLAB Online 中打开

Your code should somehow look like this. Fill in the details.

nx = 11;
ny = 11;
Lx = 0.1;
Ly = 0.1;
dx = Lx/(nx-1);
dy = Ly/(ny-1);
for ix = 1:nx
  for iy = 1:ny
    T(ix,iy) = 600;
  end
end
for ix = 1:nx
  for iy = 1:ny
    Y0((ix-1)*ny + iy ) = T(ix,iy);
  end
end
[time,Y] = ode15s(@(time,Y)fun(time,Y,nx,ny,dx,dy),tspan,Y0)
function dYdt = fun(time,Y,nx,ny,dx,dy)
  % Write solution vector Y in matrix
  T = zeros(nx,ny);
  for ix = 1:nx
    for iy = 1:ny
        T(ix,iy) = Y((ix-1)*ny + iy);
    end
  end
  % Allocate workspace for the time derivatives in the grid points
  dTdt = zeros(nx,ny);
  % Set the dTdt expressions of your attached paper (Don't use function
  % handles, but just the variables here !)
  ...
  % Write back the dTdt matrix into a column vector to return it to ODE15S
  for ix = 1:nx
    for iy = 1:ny
        dYdt((ix-1)*ny + iy) = dTdt(ix,iy);
    end
  end
end