Azzera filtri
Azzera filtri

ODE45 system of equations inaccurate?

4 visualizzazioni (ultimi 30 giorni)
So I have a system of equations that I have to solve and plot for 3 separate concentrations C1,C2,C3. The initial conditions are c1(0) = .063, c2(0) = 700, and c3(0) = 5.6
There is no way that the graphs I'm getting with this are right. I know the equations are all right but did I miss something or do something incorrectly? I attached an image of the equations.
a = 1;
b = (2*(10^-3));
d = (4*(10^-4));
f = (4*(10^-8));
c1(1) = .063;
c2(1) = 700;
c3(1) = 5.6;
c1prime = @(x1,y1) (-a*(c1^3)) + (b*(c3)) + (d*(c3^4)) - f*(c2^4)*(c1^4);
c2prime = @(x2,y2) (d*(c3^4)) - f*(c2^4)*(c1^4);
c3prime = @(x3,y3) (a*(c1^3)) - (b*c3) - (d*(c3^4)) + f*(c2^4)*(c1^4);
options = odeset('RelTol',10^-7,'AbsTol',10^-9);
[C1 Y1] = ode45(c1prime,[0 5],.063,options);
[C2 Y2] = ode45(c2prime,[0 5],700,options);
[C3 Y3] = ode45(c3prime,[0 5],5.6,options);
plot(C1,Y1,'r',C2,Y2,'b',C3,Y3,'g')

Risposta accettata

Star Strider
Star Strider il 12 Mar 2016
The error you made was in not putting them all in the same system. You did code them correctly otherwise, so all I needed to do was to copy your code to the ‘cprime’ matrix I created, subscript the concentrations, and provide the correct initial conditions vector. You have to plot them using the subplot function or you won’t be able to see the changes.
This works:
a = 1;
b = (2E-3);
d = (4E-4);
f = (4E-8);
C10 = .063;
C20 = 700;
C30 = 5.6;
cprime = @(t,C) [(-a*(C(1).^3)) + (b*(C(3))) + (d*(C(3).^4)) - f*(C(2).^4).*(C(1).^4); (d*(C(3).^4)) - f*(C(2).^4).*(C(1).^4); (a*(C(1).^3)) - (b*C(3)) - (d*(C(3).^4)) + f*(C(2).^4).*(C(1).^4)]
options = odeset('RelTol',10^-7,'AbsTol',10^-9);
[X,C] = ode45(cprime, [0 1], [C10; C20; C30], options);
figure(1)
subplot(3,1,1)
plot(X, C(:,1))
ylabel('\itC_1\rm')
grid
subplot(3,1,2)
plot(X, C(:,2))
ylabel('\itC_2\rm')
grid
subplot(3,1,3)
plot(X, C(:,3))
grid
xlabel('\itx\rm')
ylabel('\itC_3\rm')
I experimented with ode15s as well. It gives the same results as ode45.

Più risposte (0)

Tag

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by