

Hello,
I'm trying to simulate an electrical circuit whose timedomain differential equation is
V" + k1 V' +k2 V = a1 I' + a2 I
Vs(t)  Rs I = V
where V and I are the device voltage and current, respectively, and the rest are constants. Vs is the timevarying voltage source, and Rs its internal resistance. Starting from equilibrium, all values and their derivatives are set to zero.
I don't see how to formulate the problem, in particular adding the source equation into the picture.
LSODE: let's say x1 = V, x2 = I, x3 = V'; so we have
x1' = x3
x2' = ??? => I' depends on Vs'(t) ???
x3' = a1 * ??? + a2 x2  k1 x3  k2 x1
How to fill the "???"
DASSL: same variables
res1 = x1'x3
res2 = ...
res3 = x3'  a1 x2'  a2 x2  k1 x3  k2 x1
In all cases, I need a relationship around I', and would like Vs(t) to be a step source, in this case Vs'(t) would become infinite ???
Any hints ?
Regards
Pascal


On May 16, 2011, at 9:35 AM, CdeMills wrote:
> Hello,
>
> I'm trying to simulate an electrical circuit whose timedomain differential
> equation is
> V" + k1 V' +k2 V = a1 I' + a2 I
> Vs(t)  Rs I = V
> where V and I are the device voltage and current, respectively, and the rest
> are constants. Vs is the timevarying voltage source, and Rs its internal
> resistance. Starting from equilibrium, all values and their derivatives are
> set to zero.
>
> I don't see how to formulate the problem, in particular adding the source
> equation into the picture.
>
> LSODE: let's say x1 = V, x2 = I, x3 = V'; so we have
> x1' = x3
> x2' = ??? => I' depends on Vs'(t) ???
> x3' = a1 * ??? + a2 x2  k1 x3  k2 x1
> How to fill the "???"
>
> DASSL: same variables
> res1 = x1'x3
> res2 = ...
> res3 = x3'  a1 x2'  a2 x2  k1 x3  k2 x1
> In all cases, I need a relationship around I', and would like Vs(t) to be a
> step source, in this case Vs'(t) would become infinite ???
>
> Any hints ?
I would substitute I = (VsV)/Rs into your first equation and use lsode.
You would then have 2 d.e.'s:
x'(1) = x(2)
x'(2) = aVs  k1 x(2)  b x(1)
where:
a = a1/Rs and b = k2+a1/Rs
The initial conditions would be:
x(1)=0, x(2)=0
Don't worry about the infinite derivative, the imbalance in the initial voltage should take care of that.
_______________________________________________
Helpoctave mailing list
[hidden email]
https://mailman.cae.wisc.edu/listinfo/helpoctave


On 05/16/2011 04:52 PM, Liam Groener wrote:
> On May 16, 2011, at 9:35 AM, CdeMills wrote:
>
>> Hello,
>>
>> I'm trying to simulate an electrical circuit whose timedomain differential
>> equation is
>> V" + k1 V' +k2 V = a1 I' + a2 I
>> Vs(t)  Rs I = V
>> where V and I are the device voltage and current, respectively, and the rest
>> are constants. Vs is the timevarying voltage source, and Rs its internal
>> resistance. Starting from equilibrium, all values and their derivatives are
>> set to zero.
>>
>> I don't see how to formulate the problem, in particular adding the source
>> equation into the picture.
>>
>> LSODE: let's say x1 = V, x2 = I, x3 = V'; so we have
>> x1' = x3
>> x2' = ??? => I' depends on Vs'(t) ???
>> x3' = a1 * ??? + a2 x2  k1 x3  k2 x1
>> How to fill the "???"
>>
>> DASSL: same variables
>> res1 = x1'x3
>> res2 = ...
>> res3 = x3'  a1 x2'  a2 x2  k1 x3  k2 x1
>> In all cases, I need a relationship around I', and would like Vs(t) to be a
>> step source, in this case Vs'(t) would become infinite ???
>>
>> Any hints ?
>
> I would substitute I = (VsV)/Rs into your first equation and use lsode.
> You would then have 2 d.e.'s:
> x'(1) = x(2)
> x'(2) = aVs  k1 x(2)  b x(1)
> where:
> a = a1/Rs and b = k2+a1/Rs
>
> The initial conditions would be:
> x(1)=0, x(2)=0
> Don't worry about the infinite derivative, the imbalance in the initial voltage should take care of that.
>
> _______________________________________________
> Helpoctave mailing list
> [hidden email]
> https://mailman.cae.wisc.edu/listinfo/helpoctaveInteresting. This works. In any case, the three variable cannot
possibly work, since lsode is based on methods (ABM, BDF) that assume
continuous solutions, and this is all they will generate. However, the
formula Vs = V+Rs*I guarantees that one of V or I cannot possibly be
continuous at any jump discontinuity of Vs. Thus the 3variable
approach cannot possibly work in any formulation since it would mandate
that V, V' and I be continuous initially and at every interior point.
In that regard, the conditions that I and I' vanish initially make the
system inconsistent if they are meant to be used as part of an initial
value problem. Your approach ignore those constraints and banished
the potential discontinuities due to the step source on the right hand
side of the first order ode as simple jump discontinuities, so the
problem remains well posed. Here's a script that illustrates these
points with an impulse function.
Thomas Shores
% DEtest.m
clear
global k = [1,2] # [k1,k2]
global a = [2,2] # [a1,a2]
global Rs = 0.5
function retval = heaviside(t)
retval = (t==0)*0.5 + (t>0);
end
function retval = Vs(t)
retval = heaviside(t0.5)heaviside(t1);
end
function retval = Vsp(t)
retval = zeros(size(t)); # ignore any dirac delta function
end
function retval = fcn(x,t)
# here x = [V; Vp]
global k;
global a;
global Rs;
retval = [0;0];
retval(1) = x(2);
retval(2) = a(1)*(Vsp(t)x(2))/Rs + a(2)*(Vs(t)x(1))/Rs  k(1)*x(2) 
k(2)*x(1);
end
% main: give it a spin
x0 = [0;0];
t = linspace(0,4,800);
x = lsode(@fcn,x0,t);
clf
plot(t',x(:,1),"1;V(t);")
hold on, grid
plot(t',x(:,2),"2;V'(t);")
plot(t',(Vs(t)'x(:,1))/Rs,"3;I(t);")
plot(t,Vs(t),"4;Vs(t);")
_______________________________________________
Helpoctave mailing list
[hidden email]
https://mailman.cae.wisc.edu/listinfo/helpoctave


Sent from my iPad
On May 16, 2011, at 8:41 PM, Thomas Shores < [hidden email]> wrote:
> On 05/16/2011 04:52 PM, Liam Groener wrote:
>> On May 16, 2011, at 9:35 AM, CdeMills wrote:
>>
>>> Hello,
>>>
>>> I'm trying to simulate an electrical circuit whose timedomain differential
>>> equation is
>>> V" + k1 V' +k2 V = a1 I' + a2 I
>>> Vs(t)  Rs I = V
>>> where V and I are the device voltage and current, respectively, and the rest
>>> are constants. Vs is the timevarying voltage source, and Rs its internal
>>> resistance. Starting from equilibrium, all values and their derivatives are
>>> set to zero.
>>>
>>> I don't see how to formulate the problem, in particular adding the source
>>> equation into the picture.
>>>
>>> LSODE: let's say x1 = V, x2 = I, x3 = V'; so we have
>>> x1' = x3
>>> x2' = ??? => I' depends on Vs'(t) ???
>>> x3' = a1 * ??? + a2 x2  k1 x3  k2 x1
>>> How to fill the "???"
>>>
>>> DASSL: same variables
>>> res1 = x1'x3
>>> res2 = ...
>>> res3 = x3'  a1 x2'  a2 x2  k1 x3  k2 x1
>>> In all cases, I need a relationship around I', and would like Vs(t) to be a
>>> step source, in this case Vs'(t) would become infinite ???
>>>
>>> Any hints ?
>>
>> I would substitute I = (VsV)/Rs into your first equation and use lsode.
>> You would then have 2 d.e.'s:
>> x'(1) = x(2)
>> x'(2) = aVs  k1 x(2)  b x(1)
>> where:
>> a = a1/Rs and b = k2+a1/Rs
>>
>> The initial conditions would be:
>> x(1)=0, x(2)=0
>> Don't worry about the infinite derivative, the imbalance in the initial voltage should take care of that.
>>
>> _______________________________________________
>> Helpoctave mailing list
>> [hidden email]
>> https://mailman.cae.wisc.edu/listinfo/helpoctave> Interesting. This works. In any case, the three variable cannot possibly work, since lsode is based on methods (ABM, BDF) that assume continuous solutions, and this is all they will generate. However, the formula Vs = V+Rs*I guarantees that one of V or I cannot possibly be continuous at any jump discontinuity of Vs. Thus the 3variable approach cannot possibly work in any formulation since it would mandate that V, V' and I be continuous initially and at every interior point. In that regard, the conditions that I and I' vanish initially make the system inconsistent if they are meant to be used as part of an initial value problem. Your approach ignore those constraints and banished the potential discontinuities due to the step source on the right hand side of the first order ode as simple jump discontinuities, so the problem remains well posed. Here's a script that illustrates these points with an impulse function.
>
> Thomas Shores
>
> % DEtest.m
>
> clear
> global k = [1,2] # [k1,k2]
> global a = [2,2] # [a1,a2]
> global Rs = 0.5
>
> function retval = heaviside(t)
> retval = (t==0)*0.5 + (t>0);
> end
>
> function retval = Vs(t)
> retval = heaviside(t0.5)heaviside(t1);
> end
>
> function retval = Vsp(t)
> retval = zeros(size(t)); # ignore any dirac delta function
> end
>
> function retval = fcn(x,t)
> # here x = [V; Vp]
> global k;
> global a;
> global Rs;
>
> retval = [0;0];
> retval(1) = x(2);
> retval(2) = a(1)*(Vsp(t)x(2))/Rs + a(2)*(Vs(t)x(1))/Rs  k(1)*x(2)  k(2)*x(1);
> end
>
> % main: give it a spin
> x0 = [0;0];
> t = linspace(0,4,800);
> x = lsode(@fcn,x0,t);
> clf
> plot(t',x(:,1),"1;V(t);")
> hold on, grid
> plot(t',x(:,2),"2;V'(t);")
> plot(t',(Vs(t)'x(:,1))/Rs,"3;I(t);")
> plot(t,Vs(t),"4;Vs(t);")
>
Actually, I see that I forgot the I' part of the original d.e., so back to the drawing boards.
_______________________________________________
Helpoctave mailing list
[hidden email]
https://mailman.cae.wisc.edu/listinfo/helpoctave


Must Vs be a perfect step source, or could you approximate it with a very
steep slope? Presumably in reality it would take some finite time to go from
one value to the other?
You could then specify the derivative algebraically in your function
something along the lines of:
if 0 <= t < t1 Vs' = 0;
elseif t1 <= t < t2 Vs' = 1000; (or perhaps make Vs some
quadratic function during this period?)
elseif t2 <= t < tend Vs' = 0;
Just a suggestion

View this message in context: http://octave.1599824.n4.nabble.com/differentialequationtp3526715p3528765.htmlSent from the Octave  General mailing list archive at Nabble.com.
_______________________________________________
Helpoctave mailing list
[hidden email]
https://mailman.cae.wisc.edu/listinfo/helpoctave


Oops, forgot to reply to the list.
 Forwarded message 
From: Joshua Stults < [hidden email]>
Date: Tue, May 17, 2011 at 9:43 AM
Subject: Re: differential equation
To: Liam Groener < [hidden email]>
On Tue, May 17, 2011 at 1:04 AM, Liam Groener < [hidden email]> wrote:
> Actually, I see that I forgot the I' part of the original d.e., so back to the drawing boards.
Extending your original approach, you can take the derivative of the
Vs  Rs I = V equation and eliminate I' as well as I; then you have a
single second order ode in V with a forcing term that includes Vs and
Vs'. This doesn't solve the op's problem with infinite Vs' though.

Joshua Stults
Website: variousconsequences.com
Hackerspace: daytondiode.org

Joshua Stults
Website: variousconsequences.com
Hackerspace: daytondiode.org
_______________________________________________
Helpoctave mailing list
[hidden email]
https://mailman.cae.wisc.edu/listinfo/helpoctave


On 05/17/2011 05:38 AM, Richard Crozier wrote:
Must Vs be a perfect step source, or could you approximate it with a very
steep slope? Presumably in reality it would take some finite time to go from
one value to the other?
You could then specify the derivative algebraically in your function
something along the lines of:
if 0 <= t < t1 Vs' = 0;
elseif t1 <= t < t2 Vs' = 1000; (or perhaps make Vs some
quadratic function during this period?)
elseif t2 <= t < tend Vs' = 0;
Just a suggestion

View this message in context: http://octave.1599824.n4.nabble.com/differentialequationtp3526715p3528765.html
Sent from the Octave  General mailing list archive at Nabble.com.
_______________________________________________
Helpoctave mailing list
[hidden email]
https://mailman.cae.wisc.edu/listinfo/helpoctave
I have responded to this topic twice two days ago and my posts have
not yet appeared on helpoctave, so let me try one more time. This
was in response to Liam Groener's response to a prior post:
Actually, I see that I forgot the I' part of the original d.e., so back to the drawing boards.
Sorry if this is a repost, but my reply hasn't appeared on
octavehelp, so here it is (maybe again).
Drats! You're not alone. I neglected the effect of the derivatives
of Vs(t) as well, which is incorrect So it's back to applied math
101, where the real problem with this model lies. lsode handles
jump discontinuities in the rhs just fine, and produces a continuous
solution, as it must. However, it cannot handle delta functions on
the rhs, because the solution in this case *must*
also suffer a discontinuity. So the right thing to do is to restart
lsode at every jump discontinuity of Vs(t) and account for the delta
function by adding the jump to the new initial condition for the
corresponding lhs, which in this case is V'. So: Your approach is
correct in ignoring the initial conditions on I, and if you make the
modifications I suggest it will produce a correct solution. Where
Vs(t) suffers a jump discontinuity, both I(t) and V'(t) will do the
same, but V(t) remains continuous across the entire interval. I
modified my original code to yield a correct (I think) solution to
the problem. Here is the script for a sample discontinuous impulse
source term Vs. The graph is busy, so rem out the plots you don't
want to see:
% DEtest.m
clear
global k = [1,2] # [k1,k2]
global a = [2,2] # [a1,a2]
global Rs = 0.5
function retval = heaviside(t)
retval = (t==0)*0.5 + (t>0);
end
function retval = Vs(t)
retval = heaviside(t0.5)1.5*heaviside(t1)+0.5*heaviside(t2);
end
function retval = Vsp(t)
retval = zeros(size(t)); % only valid off discontinuities
end
function retval = fcn(x,t)
# here x = [V; Vp]
global k;
global a;
global Rs;
retval = [x(2);
a(1)*(Vsp(t)x(2))/Rs+a(2)*(Vs(t)x(1))/Rsk(1)*x(2)k(2)*x(1)];
end
% main: give it a spin
bdy = [0.5,1; 1,1.5; 2,0.5; 4,0]; # [point, jump] pairs for lsode
up to right bdy
xsoln = x0 = [0,0];
tsoln = 0;
for j = 1:rows(bdy)
t = linspace(tsoln(end),bdy(j,1),100);
x = lsode(@fcn,x0,t);
tsoln = [tsoln,t(2:end)];
xsoln = [xsoln;x(2:end,:)];
x0 = (xsoln(end,:)+[0,bdy(j,2)])';
end
clf
plot(tsoln',xsoln(:,1),"1;V(t);")
hold on, grid
plot(tsoln',xsoln(:,2),"2;V'(t);")
plot(tsoln',(Vs(tsoln)'xsoln(:,1))/Rs,"3;I(t);")
plot(tsoln,Vs(tsoln),"4;Vs(t);")
_______________________________________________
Helpoctave mailing list
[hidden email]
https://mailman.cae.wisc.edu/listinfo/helpoctave

