How I can to resolve the code below?

%Function
function di=didt(t,i,flag,L,R,v,f,Rsec) 
i1=i(1);
i2=i(2);
A=L;
D=L*R;
w=2*pi*f;
di1=A(1,1)*(v*cos(w*t))+A(1,2)*Rsec*i2+A(1,3)*((D(3,1)*i1+D(3,2)*i2-A(3,1)*(v*cos(w*t))-A(3,2)*Rsec*i2)/A(3,3))-D(1,1)*i1-D(1,2)*i2;
di2=A(2,1)*(v*cos(w*t))+A(2,2)*Rsec*i2+A(2,3)*((D(3,1)*i1+D(3,2)*i2-A(3,1)*(v*cos(w*t))-A(3,2)*Rsec*i2)/A(3,3))-D(2,1)*i1-D(2,2)*i2;
di=[di1 di2]';
end

_________________________________________________________________________________

%Code
R=[2.1581 0 0; 0 0.0114 0; 0 0 0.0556];
L=[2.9619 -35.9491 0.1375; -35.9491 609.1265 -174.4702; 0.1375 -174.4702 172.8074];
v=10e3;
f=1e3;
Rsec=100;
i10=0;
i20=0;
ts=0:0.01:10e-3;  
[t,i]=ode45('didt',ts,[i10 i20],[],L,R,v,f,Rsec)
figure(1)
plot(t,i(:,1), t,i(:,2))

[Image Analyst]

Resolve in what way? What's wrong with it. Please read this link.

[James Tursa]

What specific problems are you having with this? An error? Results not as expected?

[Aline Silva]

The problem of my program are results not as expected. I need to find the primary and secundary current of the transformer using the function ode45. But I don't properly implement this function and to find expected results.

[James Tursa]

What equations are you implementing? We can see a bunch of constants and some code, but without knowing what you are trying to implement it is very difficult to help you.

[Aline Silva]

I need to implement the di1 and di2 equations of the function called didt.

[James Tursa]

No, that's the code for the equations. What are the original equations? From a book? From a website? From a homework handout?

[Aline Silva]

The differential equation is: [di(t)/dt]=[L][v(t)]-[L][R][i(t)]

[Aline Silva]

The image below shows all conditions:

[Walter Roberson]

Your v_1 and v_2 appear to be formulas rather than initial conditions.

You do not appear to define v_3

Your definition of v_2 = R*i_2 = 100*i_2 would seem to imply that v_2 should be 3 x 3 since you have defined R as a 3 x 3 matrix -- but another possibility is you are defining R there as being equal to 100, in contradiction to the structural matching which would suggest that R is a 3 x 3 diagonal matrix.

Is i_3 = 0 defining i_3 as being constant even though it also appears to be a differential equation? Or should that be i_3(0) = 0 for the initial condition?

It is not clear to me where you define initial conditions for i_1 or i_2 ?

What is the value of omega for the definition of v_1 ?

[Aline Silva]

It's a transformer of three windings. The conditions initials are:

- The primary voltage is cosinusoidal with amplitude of 10 kV and frequency of 1 kHz.

- The secundary there was a resistive load of 100 Ohm: v_2=R*i_2(t) => v_2=100*i_2(t)

- The terciary is open => i_3(t)=0 (as constant).

I need to find the primary [i_1(t)] and secundary [i_2(t)] current of the transformer, solving the differential equations system.

[Walter Roberson]

An initial condition would be something like i_2(0) = 1/3, or di_1(t)/dt = -17 . Initial conditions are initial state, not equations. For example you might have the equations for a pendulum, but how the pendulum acts is going to depend upon how high you raise it at your start time. The indefinite integration of a formula has a constant, and in order to model the system properly you need the information that allows you to figure out what the proper constant is.

I still do not know what v_3 should be set to in your system.

Frequency of 1 kHz implies omega = 2*Pi*1000, so at least we know that.

Answer 1

[Walter Roberson]

   The differential equation is: [di(t)/dt]=[L][v(t)]-[L][R][i(t)]

What is v(t) here? Is it a scalar function? Is it a function that returns a vector? Is it a function that returns an array?

Your L appears to be a 3 x 3 array, and your R appears to be a 3 x 3 array as well. In order for the dimensions to work out, it would appear that you need v(t) and i(t) to be both 3 x N for some positive integer N.

[L][v(t)] -> (3 x 3) by (3 x N) -> 3 x N
[L][R][i(t)] -> (3 x 3) by (3 x 3) by (3 x N) -> 3 x N
(3 x N) minus 3 x N matches dimensions properly giving 3 x N for right side
derivative of 3 x N is 3 x N so left side is 3 x N
left side and right side match

However, we cannot tell from this analysis what N is .

You pass in a constant scalar for v to didt . Does that work out?

[di(t)/dt]=[L][v(t)]-[L][R][i(t)]

(3 x 3) by (1 x 1) -> 3 x 3  (permitted as an exception to normal matrix multiplication)
[L][R][i(t)] -> (3 x 3) by (3 x 3) by Something -> (3 x 3) by Something. This can only work if Something is a scalar, which gives 3 x 3, or is 3 x N, giving 3 x N
In the case that the Something is a scalar, (3 x 3) minus (3 x 3) -> 3 x 3 works okay giving 3 x 3. So hypothetically, i(t) might be a scalar up to this point, with the right side yielding 3 x 3
In the case that the Something is 3 x N, (3 x 3) minus (3 x N). Before R2016b this only works okay if N is 3, in which case it would be (3 x 3) minus (3 x 3) giving 3 x 3. In R2016b, (3 x 3) minus (3 x N) still works okay if N is 3, giving (3 x 3), but in R2016b, (3 x 3) minus (3 x 1) also works, also giving (3 x 3). So hypothetically, i(t) might be 3 x 3 (valid in all versions) or might be 3 x 1 (only in R2016b or later); either way the right side would yield 3 x 3
derivative of i(t) is derivative of size Something. Under the hypothesis that i(t) was scalar, di(t)/dt would be a scalar as well so the left would be scalar under the hypothesis that i(t) is scalar. Under the hypothesis, valid in versions so far, that i(t) was 3 x 3, di(t)/dt would be 3 x 3. Under the hypothesis valid in R2016b and later only, that i(t) was 3 x 1, di(t)/dt would be 3 x 1.
The size of the left must equal the size of the right. Under all of the hypotheses so far abotu the size of i(t), the right side must come out as 3 x 3. Under the hypothesis that i(t) was a scalar, the left was a scalar and so does not match the 3 x 3 right, so we must reject the hypothesis that i(t) is a scalar. Under the hypothesis, valid only for R2016b and later, that i(t) was 3 x 1, the left was 3 x 1 and so does not match the 3 x 3 right, so we must reject the hypothesis that i(t) is 3 x 1. Under the hypothesis, valid for all releases, that i(t) was 3 x 3, the left was 3 x 3 and that _does_ match the 3 x 3 right. Having exhausted all other hypotheses, we must conclude that if v(t) is a scalar then i(t) is 3 x 3.

We now have two competing hypothesis: that v(t) is a scalar, in which case i(t) must be 3 x 3; or that v(t) and i(t) are both 3 x N for some as-yet undetermined positive integer N

Under either hypotheses, the number of rows in i(t) must be 3, so the minimum number of items that di(t)/dt could return for an ode would have to be 3. There must be at least one original condition for each item in i(t), the boundary conditions; there could also be additional boundary conditions having to do with second derivatives and so on. Your code passes in two boundary conditions total, which cannot be enough to satisfy the i(t) that we have proven to require at least 3 under the assumption that those L and R in the code are the correct size.

At this point, we are driven to one of three hypotheses:

1. the L and R given in your code are the wrong size and the correct sizes are 1 by something or 2 by something; and additionally that your didt code is wrong in referring to those non-existent values; this hypothesis would allow the initial conditions to match; Or,
2. your didt code is wrong in not accepting at least 3 or at least 9 initial conditions and returning the same number of values; if the i(t) is not 3 x 3 then, additionally, that the size of v is wrong and must be 3 x N for some positive integer N; this hypothesis would allow the initial conditions to match without changing the size of L or R; Or,
3. your given differential equation does not correctly express your hypothetically correct code

[Aline Silva]

The image below shows all conditions.

[Walter Roberson]

Your v_1 and v_2 appear to be formula rather than initial conditions.

You do not define v_3

Your definition of v_2 implies that v_2 is 3 x 3 since R is 3 x 3 and that is multiplied by the scalar i_2. But then you have the = 100i_2 which implies that maybe R = 100 so that Ri_2 and 100i_2 can be the same. But if so then what does that mean about the R in the rest of your equations?

Are you defining i_3(t) as the constant 0? Or are you defining i_3(0) as the initial condition 0? What are the initial conditions for i_1(t) and i_2(t) ?

[Aline Silva]

It's a transformer of three windings. The conditions initials are:

- The primary voltage is cosinusoidal with amplitude of 10 kV and frequency of 1 kHz.

- The secundary there was a resistive load of 100 Ohm: v_2=R*i_2(t) => v_2=100*i_2(t)

- The terciary is open => i_3(t)=0 (as constant).

I need to find the primary [i_1(t)] and secundary [i_2(t)] current of the transformer, solving the differential equations system.

[Walter Roberson]

If i_3(t) is 0 then we can use the fact that the derivative of i_3(t) would have to be 0 as well, and put that 0 on the left side of the equation for di_3(t)/dt . That gives us a non-differential equation that we can solve for v_3, getting that

 v_3 = -(2712875/340949)*cos(omega*t) + (860476156034189/8523725000000)*i_2(t) + (66910349/38965600000)*i_1(t)

We can then substitute that into the first two differential equations, along with i_3(t) = 0, to remove the v_3 from them. Substitute in the omega, and let i_1(0) = init_i1 (a symbolic constant) and let i_2(0) = init_i2 (a symbolic constant) and Maple's dsolve() can solve, getting

i_1(t) = ((((1054921529115435732580778844485540706364492187500000000000000000000000000000000000000000000*init_i1+174473063835960627181610309086312603007812500000000000000000000000000000000000000000000000*init_i2)*Pi^4+(494155401516438578715044412561608046862245124595971314286475475502192044326171875000000000000*init_i1+81728170801439764518912993617767647661926523934585016131717904140582548828125000000000000000*init_i2-338013189172899933695684825820411936260700409632954301129966903686523437500000000000000000000)*Pi^2+499128614813175122702845533033732067862732520236968875754094847276124629212739*init_i1+82550688625792949458623085563182725542319825973928564408758314313895075453000*init_i2-2312815044776308431967219002982864871241983783128533783207890492915641671900000000)*3737875373307147684900657499628282286652828915609^(1/2)+2038939433217361816816712552605626392546248446773069326970946411835937500000000000000000000000000000000000000000000*Pi^4*init_i1+(955097518138690626756144458987301741585560956456426804229742936629484806388687939520704146703185224609375000000000000*init_i1+653500697206873893954480331148629818897221894281453433381966601258704587414367484019973754882812500000000000000000000)*Pi^2+964709683992410225097425816947322181949615807028424031369590093925721971656867963301142326837667138367*init_i1-4470180640342941592592677896980316861821119535834410042952551290142819941878819161767954806717330700000000)*exp(-(1/82093018322000000000000)*(-1776756523728883013591422157+3156863763155657955907384198543561716297799841795652649^(1/2))*t)+(((-1054921529115435732580778844485540706364492187500000000000000000000000000000000000000000000*init_i1-174473063835960627181610309086312603007812500000000000000000000000000000000000000000000000*init_i2)*Pi^4+(-494155401516438578715044412561608046862245124595971314286475475502192044326171875000000000000*init_i1-81728170801439764518912993617767647661926523934585016131717904140582548828125000000000000000*init_i2+338013189172899933695684825820411936260700409632954301129966903686523437500000000000000000000)*Pi^2-499128614813175122702845533033732067862732520236968875754094847276124629212739*init_i1-82550688625792949458623085563182725542319825973928564408758314313895075453000*init_i2+2312815044776308431967219002982864871241983783128533783207890492915641671900000000)*3737875373307147684900657499628282286652828915609^(1/2)+2038939433217361816816712552605626392546248446773069326970946411835937500000000000000000000000000000000000000000000*Pi^4*init_i1+(955097518138690626756144458987301741585560956456426804229742936629484806388687939520704146703185224609375000000000000*init_i1+653500697206873893954480331148629818897221894281453433381966601258704587414367484019973754882812500000000000000000000)*Pi^2+964709683992410225097425816947322181949615807028424031369590093925721971656867963301142326837667138367*init_i1-4470180640342941592592677896980316861821119535834410042952551290142819941878819161767954806717330700000000)*exp((1/82093018322000000000000)*(1776756523728883013591422157+3156863763155657955907384198543561716297799841795652649^(1/2))*t)+(-1307001394413747787908960662297259637794443788562906866763933202517409174828734968039947509765625000000000000000000000*Pi^2+8940361280685883185185355793960633723642239071668820085905102580285639883757638323535909613434661400000000)*cos(2000*Pi*t)+60389116340524053425071561315836121699438192875764362987972645239171135253906250000000000000000000000000000000000000*sin(2000*Pi*t)*Pi*(Pi^2+15053930646795331509620433548259289061977550931/3231200096812953601871582031250000000000000000000))/(4077878866434723633633425105211252785092496893546138653941892823671875000000000000000000000000000000000000000000000*Pi^4+1910195036277381253512288917974603483171121912912853608459485873258969612777375879041408293406370449218750000000000000*Pi^2+1929419367984820450194851633894644363899231614056848062739180187851443943313735926602284653675334276734)

i_2(t) = (1054921529115435732580778844485540706364492187500000000000000000000000000000000000000000000*(3156863763155657955907384198543561716297799841795652649^(1/2)-1777281249190089455846595907)*(((init_i1+(319851859252535019531000/1933929542100206154348853)*init_i2)*Pi^4+((3737875362329543267593204395324100470626643955609/7979605566935044000000000000000000000000000000)*init_i1+(9126460948819445896114166353827052476097356653109602673814518412209/117801488093152256288191661587044022172000000000000000000000000000)*init_i2-184405107641372631087768784244728634318523413102011147193/575518720455094992552997430654260871680000000000000000)*Pi^2+(36870054610686009456688014330771178647941100898785009/77925835614600039062500000000000000000000000000000000000000000000)*init_i1+(90022561205881079017037170734114204517251718619333221819801869480801609/1150405157159690002814371695185976779023437500000000000000000000000000000000000000000)*init_i2-11959148430323065020009086711245922363912131332723/5454808493022002734375000000000000000000000000000000000000)*3737875373307147684900657499628282286652828915609^(1/2)+(3737875373307147684900657499628282286652828915609/1933929542100206154348853)*Pi^4*init_i1+((13971712265343131454261313228614242619765920277305443536160810958694177739171663572417157803200881/15431994940202945573753107667902766904532000000000000000000000000000000)*init_i1+356521195835654661434988162379453443532664047746379516762431120231841832733590179/575518720455094992552997430654260871680000000000000000)*Pi^2+(137815669141772889299632259563903168849945115289774861624227156275103138808123994757306046691095305481/150703075587919390368682692069362958052070312500000000000000000000000000000000000000000000)*init_i1-44701806403429415925926778969803168618211195358344100429525512901428199418788191617679548067173307/10549215291154357325807788444855407063644921875000000000000000000000000000000000000)*exp(-(1/82093018322000000000000)*(-1776756523728883013591422157+3156863763155657955907384198543561716297799841795652649^(1/2))*t)+1054921529115435732580778844485540706364492187500000000000000000000000000000000000000000000*(3156863763155657955907384198543561716297799841795652649^(1/2)+1777281249190089455846595907)*(((init_i1+(319851859252535019531000/1933929542100206154348853)*init_i2)*Pi^4+((3737875362329543267593204395324100470626643955609/7979605566935044000000000000000000000000000000)*init_i1+(9126460948819445896114166353827052476097356653109602673814518412209/117801488093152256288191661587044022172000000000000000000000000000)*init_i2-184405107641372631087768784244728634318523413102011147193/575518720455094992552997430654260871680000000000000000)*Pi^2+(36870054610686009456688014330771178647941100898785009/77925835614600039062500000000000000000000000000000000000000000000)*init_i1+(90022561205881079017037170734114204517251718619333221819801869480801609/1150405157159690002814371695185976779023437500000000000000000000000000000000000000000)*init_i2-11959148430323065020009086711245922363912131332723/5454808493022002734375000000000000000000000000000000000000)*3737875373307147684900657499628282286652828915609^(1/2)-(3737875373307147684900657499628282286652828915609/1933929542100206154348853)*Pi^4*init_i1+(-(13971712265343131454261313228614242619765920277305443536160810958694177739171663572417157803200881/15431994940202945573753107667902766904532000000000000000000000000000000)*init_i1-356521195835654661434988162379453443532664047746379516762431120231841832733590179/575518720455094992552997430654260871680000000000000000)*Pi^2-(137815669141772889299632259563903168849945115289774861624227156275103138808123994757306046691095305481/150703075587919390368682692069362958052070312500000000000000000000000000000000000000000000)*init_i1+44701806403429415925926778969803168618211195358344100429525512901428199418788191617679548067173307/10549215291154357325807788444855407063644921875000000000000000000000000000000000000)*exp((1/82093018322000000000000)*(1776756523728883013591422157+3156863763155657955907384198543561716297799841795652649^(1/2))*t)-214623066543898322185274902482296800288196923177738070275674272885776169024407130946456054332774414062500000000000000000000000000000000000000000*((Pi^2+192015766567972249565422503/279151993750000000000000000000000)*sin(2000*Pi*t)-(1933358567713692071372603/89328638000000000000000)*Pi*cos(2000*Pi*t))*Pi)/(1198667449119669207229559505170075267303665392621111078037850503366360650421123126187298984375000000000000000000000000000000000000000000000000*Pi^4+561490099743332972457042945078503281895909324115233842819240655620770138114367178740040985255887042786174585638896484671152343750000000000000000*Pi^2+567140973985444397686736818049492647173552883747884073686341645102432195801885983687795004214275982088974551671842518531521926000)

More compactly, optimized code for this is

t11 = init_i1 + 319851859252535019531000 ./ 1933929542100206154348853 * init_i2;
t2 = 1054921529115435732580778844485540706364492187500000000000000000000000000000000000000000000 * t11;
t3 = 494155401516438578715044412561608046862245124595971314286475475502192044326171875000000000000 * init_i1 + 81728170801439764518912993617767647661926523934585016131717904140582548828125000000000000000 * init_i2 - 338013189172899933695684825820411936260700409632954301129966903686523437500000000000000000000;
t4 = pi ^ 2;
t51 = t4 ^ 2;
t6 = 82550688625792949458623085563182725542319825973928564408758314313895075453000 * init_i2;
t7 = 499128614813175122702845533033732067862732520236968875754094847276124629212739 * init_i1;
t8 = sqrt(3737875373307147684900657499628282286652828915609);
t91 = 955097518138690626756144458987301741585560956456426804229742936629484806388687939520704146703185224609375000000000000 * init_i1 + 653500697206873893954480331148629818897221894281453433381966601258704587414367484019973754882812500000000000000000000;
t10 = t91 * t4;
t11 = (2038939433217361816816712552605626392546248446773069326970946411835937500000000000000000000000000000000000000000000 * t51 + 964709683992410225097425816947322181949615807028424031369590093925721971656867963301142326837667138367) * init_i1;
t12 = sqrt(3156863763155657955907384198543561716297799841795652649);
t7 = t12 / 82093018322000000000000;
t13 = exp(-(-1933358567713692071372603 ./ 89328638000000000000 + t7) * t);
t14 = exp((1933358567713692071372603 ./ 89328638000000000000 + t7) * t);
t151 = 2000 * pi * t;
t16 = cos(t151);
t152 = sin(t151);
t17 = 60389116340524053425071561315836121699438192875764362987972645239171135253906250000000000000000000000000000000000000 * t152 * pi * (t4 + 15053930646795331509620433548259289061977550931 ./ 3231200096812953601871582031250000000000000000000);
t18 = t4 * (4077878866434723633633425105211252785092496893546138653941892823671875000000000000000000000000000000000000000000000 * t4 + 1910195036277381253512288917974603483171121912912853608459485873258969612777375879041408293406370449218750000000000000) + 1929419367984820450194851633894644363899231614056848062739180187851443943313735926602284653675334276734;
t12 = (36870054610686009456688014330771178647941100898785009 ./ 77925835614600039062500000000000000000000000000000000000000000000 * init_i1 + 90022561205881079017037170734114204517251718619333221819801869480801609 ./ 1150405157159690002814371695185976779023437500000000000000000000000000000000000000000 * init_i2 + t4 * (t11 * t4 + t3 / 1054921529115435732580778844485540706364492187500000000000000000000000000000000000000000000) - 11959148430323065020009086711245922363912131332723 ./ 5454808493022002734375000000000000000000000000000000000000) * t8;
t52 = (3737875373307147684900657499628282286652828915609 ./ 1933929542100206154348853 * t51 + 137815669141772889299632259563903168849945115289774861624227156275103138808123994757306046691095305481 ./ 150703075587919390368682692069362958052070312500000000000000000000000000000000000000000000) * init_i1;
t15 = (-214623066543898322185274902482296800288196923177738070275674272885776169024407130946456054332774414062500000000000000000000000000000000000000000 * t152 * (192015766567972249565422503 ./ 279151993750000000000000000000000 + t4) + 4645132331824332569582916940942593595482638156116167980596029347526774333713375998464855846930437515966921027029785156250000000000000000000000000 * pi * t16) * pi;
t37 = t91 * t4 / 1054921529115435732580778844485540706364492187500000000000000000000000000000000000000000000;
t1 = 1054921529115435732580778844485540706364492187500000000000000000000000000000000000000000000 * (-1777281249190089455846595907 + t12) * (t37 + t12 + t52 - 44701806403429415925926778969803168618211195358344100429525512901428199418788191617679548067173307 ./ 10549215291154357325807788444855407063644921875000000000000000000000000000000000000) * t13 + 1054921529115435732580778844485540706364492187500000000000000000000000000000000000000000000 * (1777281249190089455846595907 + t12) * (-t37 + t12 + 44701806403429415925926778969803168618211195358344100429525512901428199418788191617679548067173307 ./ 10549215291154357325807788444855407063644921875000000000000000000000000000000000000 - t52) * t14 + t15;
t9 = 1 ./ (t18);
t5 = 1 ./ (t18) / 293943858653079682948989000;
t48 = t2 * t4 + t3;
i_1(t) = ((t11 - 4470180640342941592592677896980316861821119535834410042952551290142819941878819161767954806717330700000000 + (t4 * t48 + t6 + t7 - 2312815044776308431967219002982864871241983783128533783207890492915641671900000000) * t8 + t10) * t13 + (t11 - 4470180640342941592592677896980316861821119535834410042952551290142819941878819161767954806717330700000000 + t8 * (-t4 * t48 - t6 - t7 + 2312815044776308431967219002982864871241983783128533783207890492915641671900000000) + t10) * t14 + (-1307001394413747787908960662297259637794443788562906866763933202517409174828734968039947509765625000000000000000000000 * t4 + 8940361280685883185185355793960633723642239071668820085905102580285639883757638323535909613434661400000000) * t16 + t17) * t9;
i_2(t) = t1 * t5;

We cannot plot these because we do not know the initial conditions init_i1 and init_i2

[Walter Roberson]

If we use init_i1 = 0 and init_i2 = 0, then by 1 ms the values have gone to +/- 10^18 and by 10 ms they are at +/- 10^183.

This suggests that either initial conditions of 0 are unstable or else that i_3(t) = 0 is not a correct equation.