Re: differential equation

[Thomas Shores]

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/differential-equation-tp3526715p3528765.html
Sent from the Octave - General mailing list archive at Nabble.com.
_______________________________________________
Help-octave mailing list
address@hidden
https://mailman.cae.wisc.edu/listinfo/help-octave

I have responded to this topic twice two days ago and my posts have not yet appeared on help-octave, 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 octave-help, 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(t-0.5)-1.5*heaviside(t-1)+0.5*heaviside(t-2);
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))/Rs-k(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);")