Matlab akuten

Eftersom vi alla älskar matlab så otroligt mycket så tycker jag att vi startar en separat gnälla-om-hjälp-tråd för detta vidunder.

Jag börjar gnälla:

Hur får man till den förbannade newton-filen i labb 2?! Vad ska stå i loopen etc...

Hade gärna hjälp mig men mitt program plotar inte ett skit at the moment, så jag skulle nog bara förvärra saken just nu.

Länk till uppgiften, så kan vi äldre försöka hjälpa också

Här:
http://www.math.chalmers.se/Math/Grundutb/CTH/tmv036a/0809/Laboration_2_2008.pdf

Problemet löst och chockerande nog så löste Jauquline och jag det helt själva *stolt* Om du fortfarande har problem Björn (eller någon annan) är det bara att maila mig.

Nä vi hade Alexs och Yvettes matlab lärare i dag eftersom vår vanliga var sjuk så problemet löstest jätte fort

Bjørn wrote:Nä vi hade Alexs och Yvettes matlab lärare i dag eftersom vår vanliga var sjuk så problemet löstest jätte fort

...och Junkos!
...och Jennys!

Om jag fattat det rätt kan man få tag på matlab till sin dator här: http://licenser.gu.se/student/matlab.htm

Inte matlab, men licensen att köra matlab hittar du där

Matlab kan man enklast ladda ner genom att googla på msdnaa +chalmers och sedan rimligtvis gå in på sökresultatet som hänvisar till en chalmerssida (ungefär första bästa resultatet). Där finns en länk för att ladda ner Matlab.

Matlab is the shit! =) Kul att se så många engagerade.

Jag behöver hjälp med en differentialekvation.

Min differentialekvation wrote:

%funktionsprogram
function func = lignin(t,L)

global NaOH HS A Ea R m n s Tmax Tstart

func=(-A.*(exp(-(Ea./(R.*((Tmax-(Tmax-Tstart).*exp(-t./(15.3-0.021.*Tmax)))))))).*(L.^n).*(NaOH.^m).*(1+(HS.^s)));

func = func'; %change output of lignin to column vector

Mitt huvudprogram wrote:
%Huvudprogram
clear all
clc
global NaOH HS A Ea R m n s Tmax Tstart t L

A=1.05e13;
Ea=1.3e5;
n=0.94;
m=0.38;
s=0.097;
R=8.314;
EA = 25;
L0=28;
%Tmax=170+273.15;
Tmax=443.15;
Tstart=293.15;
t=[0:20:200];
%L=4.5/28;

ConsumedNaOH = [0 0.285 .34 .385 .415 .44 .45 .45 .46 .46 .465];
NaOH = (20.59-(ConsumedNaOH.*39.9971))./(0.;
HS = (2.*(EA-NaOH.*0.).*((33./80)./0.;

[t,L] = ode15s('lignin',[0:20:200],L0);


subplot(1,1,1);
plot(t,L)
xlabel('x');
ylabel('y');

Matlab svarar med

Matlab wrote:
??? Error using ==> funfun\private\odearguments
Solving LIGNIN requires an initial condition vector of length 1.

Error in ==> ode15s at 228
[neq, tspan, ntspan, next, t0, tfinal, tdir, y0, f0, odeArgs, odeFcn, ...

Error in ==> Headprog at 24
[t,L] = ode15s('lignin',[0:20:200],L0); 11E216E0ut

Jag anser att L0 i högsta grad är en vektor med längden 1. Detta ger alltså 2 frågor

1) Kommer någon läsa denna tråd innan uppgiften ska vara inlämnad?
2) Kan denna någon isåfall hitta felet?

jag ska titta på den.

när ska den vara inne?

Vilka svar vill du få? (om du vet)

Problemet är följande: Jag vet inte riktigt hur uppgiften ser ut. MEN!

(har bytt namn på din tidsvektor, nu heter den t_span=[0:20:200], efterson du har skrivt att ODElösaren ska skriva lösningar i t...[t,L]=ode15s(....))

det största problemet är att för varje tid t_span spottar lignin(t_span,L) ut en vektor... dv
s för en skalär (t_span) kommer en vektor ut.

Code:


>> lignin(20,28)

ans =

  -0.5510
  -0.4343
  -0.3942
  -0.3539
  -0.3213
  -0.2887
  -0.2737
  -0.2737
  -0.2570
  -0.2570
  -0.2480

>> 


Stackars ODE har alltså ropat på funktionen med ett t_span och ett L och får ut en vektor... när den egentligen borde få ut ett nytt värde på L för just det t_span, som den sen ska spara i sin [t,L] vektor.

Detta beror i sin tur på att för varje tid t_span använder lignin(t_span,L) en vektor i sin beräkning, nämligen "NaOH" och "HS" (och kanske något mer?), som skapas i ode-lösnings-filen enligt följande (och görs tillgängligt genom global)

Code:

ConsumedNaOH = [0 0.285 .34 .385 .415 .44 .45 .45 .46 .46 .465];
NaOH = (20.59-(ConsumedNaOH.*39.9971))./(0.8);

alltså:

Code:


>> size(NaOH)

ans =

    1    11

>> size(lignin(20,28))

ans =

    11    1

>>

Lösningen på detta skulle kunna vara att bara ha ett värde på ConsumedNaOH (förutsättning för att använda ODE-lösare på detta viset).

om jag använder ett värde på NaOH i följande ode-lösningsfil:

Code:


%Huvudprogram
clear all
clc
global NaOH HS A Ea R m n s Tmax Tstart t L

A=1.05e13;
Ea=1.3e5;
n=0.94;
m=0.38;
s=0.097;
R=8.314;
EA = 25;
L0=28;
%Tmax=170+273.15;
Tmax=443.15;
Tstart=293.15;
t_span=[0:20:200];
%L=4.5/28;

%ConsumedNaOH = [0 0.285 .34 .385 .415 .44 .45 .45 .46 .46 .465];
ConsumedNaOH = .385;
NaOH = (20.59-(ConsumedNaOH.*39.9971))./(0.8);
HS = (2.*(EA-NaOH.*0.8)).*((33./80)./0.8);

% ODE
[t,L] = ode15s('lignin',[0:20:200],L0)


subplot(1,1,1);
plot(t,L)
xlabel('x');
ylabel('y');


får jag följande svar:

Code:



t =

    0
    20
    40
    60
    80
  100
  120
  140
  160
  180
  200


L =

  28.0000
  25.8796
  18.2697
  12.2147
    8.0645
    5.2710
    3.4114
    2.1836
    1.3796
    0.8605
    0.5294

>>


se graf här.


Alternativt, om frågan är sådan att:

Räkna ut slutkonventrationen av bla bla bla om startkoncentrationen är L0. Var tjugonde sekund togs förbrukas följande mängd NaOH (11 värden).

I så fall behövs ingen ODE lösare, endast en eller två for-loopar.

------------------------------------- EDIT ---------------------------------------------------------------------
hmm den lösingen med looparna fungerade inte. Något i unktionen ligning(t_span,L) gör att det blir imaginära tal. Orkade inte bena ut vad problemet var...

men däremot tror jag att jag har hittat lösningen...

Har kommenterat i filen vad som gör vad och varför, det borde vara hyfsat tydligt...

Code:

%Huvudprogram
clear all
clc
global NaOH HS A Ea R m n s Tmax Tstart

A=1.05e13;
Ea=1.3e5;
n=0.94;
m=0.38;
s=0.097;
R=8.314;
EA = 25;
L0=28;
%Tmax=170+273.15;
Tmax=443.15;
Tstart=293.15;
t_span=[0:20:200];
%L=4.5/28;

ConsumedNaOH_vec = [0 0.285 .34 .385 .415 .44 .45 .45 .46 .46 .465];
% NaOH = (20.59-(ConsumedNaOH.*39.9971))./(0.8);
% HS = (2.*(EA-NaOH.*0.8)).*((33./80)./0.8);


% % Lösning med for-loop
% for i=1:1:length(ConsumedNaOH_vec)
%    NaOH = (20.59-(ConsumedNaOH_vec(3).*39.9971))./(0.8);  % nytt värde på ConsumedNaOH väljs för varje i
%    HS = (2.*(EA-NaOH.*0.8)).*((33./80)./0.8);          % nytt HS beräknas för varje nytt NaOH
%    t(i)=t_span(i);                                    % varje nytt värde t_span som väljs sparas i t
%    L_save(i)=lignin(t_span(i),L0);                          % Nytt L evaluers för varje L0
%    L0=L_save(length(L_save));                                    % L0 erhåller ett nytt värde som är det sista sparade i vektorn L
% end
%
% subplot(1,1,1);
% plot(t,L_save)
% xlabel('x');
% ylabel('y');
%
% % Något skumt...  Imaginära delar i L_save


% ODE
for i=1:1:length(ConsumedNaOH_vec)                      % Kör loopen så många gånger som det finns värden i ConsumedNaOH_vec
    ConsumedNaOH=ConsumedNaOH_vec(i);                  % ConsumeNaOH får erhålla ett värde åt gången från ConsumedNaOH_vec
    NaOH = (20.59-(ConsumedNaOH.*39.9971))./(0.8);      % NaOH evalueras för det värdet ConsumedNaOH
HS = (2.*(EA-NaOH.*0.8)).*((33./80)./0.8);              % HS evalueras för värdet NaOH
[t,L] = ode15s('lignin',t_span,L0);                    % Hela tidsprocessen löses för just det ConsumedNaOH, L erhålls.
L_save(:,i)=L;                                          % Det erhållna L sparas i L_save, i en viss kolumn.
                            % Då loopen startar om väljs ett nytt värde för
                            % ConsumedNaOH, och hela tidsprocessen löses
                            % för det nya värdet. Lösningarna sparas i
                            % nästa egna vektorn i L_save.
end

L_save
t
diag(L_save)


% Eftersom vi är intresserade av lösningen där, för vatje tids-steg, väljs
% ett nytt ConsumedNaOH, plottas alltså DIAGONALEN(viktigt!!) av den
% KVADRATISKA(viktigt!) matrisen L_save, mot tiden t. Tidsvektorn har lika
% många rader som kolumner och rader i L_save.

% Vill man ha en finare graf och därmed minskar steglängden i t_span, ökar
% antalet värden i den vektorn. Det måste kontras av att ConsumedNaOH_vec
% får också lika många värden.

format compact
L_save
t
diag(L_save)

subplot(1,1,1);
plot(t,diag(L_save))
xlabel('x');
ylabel('y');




som svar får jag följande:

Code:


L_save =
  28.000  28.000  28.000  28.000  28.000  28.000  28.000  28.000  28.000  28.000  28.000
  24.801  25.427  25.649  25.880  26.059  26.241  26.332  26.332  26.415  26.415  26.463
  14.330  16.514  17.336  18.270  19.026  19.776  20.134  20.134  20.560  20.560  20.781
  7.5628  10.029  11.062  12.215  13.251  14.317  14.791  14.791  15.449  15.449  15.787
  3.8821  5.9809  6.9440  8.0645  9.1094  10.241  10.774  10.774  11.494  11.494  11.881
  1.9391  3.5144  4.3105  5.2710  6.2111  7.2796  7.8038  7.8038  8.5055  8.5055  8.8842
  0.9387  2.0313  2.6402  3.4114  4.1998  5.1395  5.6199  5.6199  6.2652  6.2652  6.6182
  0.4397  1.1514  1.5901  2.1836  2.8160  3.6051  4.0237  4.0237  4.5903  4.5903  4.9083
  0.1985  0.6398  0.9435  1.3796  1.8697  2.5103  2.8621  2.8621  3.3422  3.3422  3.6213
  0.0862  0.3479  0.5502  0.8605  1.2283  1.7338  2.0214  2.0214  2.4176  2.4176  2.6553
  0.0358  0.1849  0.3146  0.5294  0.7981  1.1869  1.4166  1.4166  1.7378  1.7378  1.9330


t =

    0
    20
    40
    60
    80
  100
  120
  140
  160
  180
  200


diag(L_save)

  28.0000
  25.4273
  17.3355
  12.2147
    9.1094
    7.2796
    5.6199
    4.0237
    3.3422
    2.4176
    1.9330

>>