codice con ode45

Sto provando a risolvere un problema di carbocementazione usando ODE45, in particolr modo sto usando il metodo delle linee. Lo scopo dell'esercizio è quello di ricavare dopo un determinato istante di tempo T il grado di assorbimento del carbonio. Non riesco proprio a capire come risolvere perchè mi sembra che sia inserito tutto correttamente. qualcuno mi può aiutare?

clear

clc

%cabocemntazione

%defnizione parametri

L=1; Temp=1273; D0=20*10^-6; Q=142000; R=8.314; T=10;

nx=100; nt=100;

%controllo parametri

if any([L Temp T D0 nx-2 nt-2])<=0

error('controlla i parametri')

end

%inizializzazione

D=D0*exp(-Q/(R*Temp)); %diffusività

dx=L/nx; %ho discretizzato dx perchè mi sefve nella formula di p

p=D/dx^2;

x=linspace(0,L,nx+1);%non ho fatto il tempo non mi serve

A=[p*ones(nx,1) -2*p*ones(nx,1) p*ones(nx,1)];

AA=spdiags(A, -1:1, nx,nx);

A(end,end-1)=2*p;

a=sparse(nx,1);

%condizione iniziale

c=GetPhi(nx);

[t,y]=ode45(@system,[0 T], c(2:end),[],A,a,p);

c(2:end)=y(end,:);

c(1)=Getg1;

plot(x,c,'b*');

axis([0 L 0 Getg1]);

xlabel('x'); ylabel('%w di C');

legend('MOL');

grid

pause(0.1)

%---------subfunction----------

function phi=GetPhi(nx)

c0=0.2;

phi=c0*ones(nx+1,1);

end

function g1=Getg1

cs=1;

g1=cs;

end

function g2=Getg2

g2=0;

end

function cdot=system(t,c,A,a,p);

a(1)=p*Getg1;

a(end)=p*Getg2;

cdot=A.*c+a ;

end

clear

clc

%cabocemntazione

%defnizione parametri

L=1; Temp=1273; D0=20*10^-6; Q=142000; R=8.314; T=10;

nx=100; nt=100;

%controllo parametri

if any([L Temp T D0 nx-2 nt-2])<=0

error('controlla i parametri')

end

%inizializzazione

D=D0*exp(-Q/(R*Temp)); %diffusività

dx=L/nx; %ho discretizzato dx perchè mi sefve nella formula di p

p=D/dx^2;

x=linspace(0,L,nx+1);%non ho fatto il tempo non mi serve

A=[p*ones(nx,1) -2*p*ones(nx,1) p*ones(nx,1)];

AA=spdiags(A, -1:1, nx,nx);

A(end,end-1)=2*p;

a=sparse(nx,1);

%condizione iniziale

c=GetPhi(nx);

[t,y]=ode45(@system,[0 T], c(2:end),[],A,a,p);

c(2:end)=y(end,:);

c(1)=Getg1;

plot(x,c,'b*');

axis([0 L 0 Getg1]);

xlabel('x'); ylabel('%w di C');

legend('MOL');

grid

pause(0.1)

%---------subfunction----------

function phi=GetPhi(nx)

c0=0.2;

phi=c0*ones(nx+1,1);

end

function g1=Getg1

cs=1;

g1=cs;

end

function g2=Getg2

g2=0;

end

function cdot=system(t,c,A,a,p);

a(1)=p*Getg1;

a(end)=p*Getg2;

cdot=A.*c+[p*Getg1; a; p*Getg2] ;

end

Answer 1

[Varun]

Hey! I ran the code you shared on my end and was able to reproduce the error. While I am unsure about what exactly you are trying to achieve, I see that the reason behind this error is that the function "system" is expected to run a column vector of dimensions [n x 1]. However, in the case the "cdot" variable is actually of size 100 x 3, which results in this error.

So, refactoring your code to ensure that the "system" method returns a column vector should resolve this error. Hope this helps!