Filter löschen
Filter löschen

ODE45 system of equations inaccurate?

4 Ansichten (letzte 30 Tage)
Kaylene Widdoes
Kaylene Widdoes am 12 Mär. 2016
Beantwortet: Star Strider am 12 Mär. 2016
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')

Akzeptierte Antwort

Star Strider
Star Strider am 12 Mär. 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.

Weitere Antworten (0)

Tags

Community Treasure Hunt

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

Start Hunting!

Translated by