How can I specify different entires in a vector - partially using a function

Hi!

I'm trying to produce the initial conditions (using a vector) for a differential equation where the first and last set of entries are 0, the middle is nonzero. The function is:

IC(x) = (x-x0)^2 (x-x1)^2 for x in (x0,x1) = 0 elsewhere

I'm thinking I can define the vector as IC(x), then redefine the rest of the vector entries part by part as zeros. In other words, if the input is [0,1] with delta x = 0.1 and x0 = 1/3, x1 = 2/3, then

IC(0) = IC(.1) = IC(.2) = IC(.3) = 0

IC(.4) = (.4-x0)^2 (.4-x1)^2 IC(.5) = (.5-x0)^2 (.5-x1)^2 IC(.6) = (.6-x0)^2 (.6-x1)^2

IC(.7) = IC(.8) = IC(.9) = IC(1) = 0

This is the code I have written so far. (The variables a & b are there so I don't get a 'dimensions don't agree' error) I can define the function, but I don't know how to modify the first and last set of entires so that they are zero. Any input will be very much appreciated!

Thank you!

E

function[IC] = Initial_Condition(L, Delta_x)

x = 0:Delta_x:L; n = length(x); x0 = floor(n/3)*Delta_x; x1 = floor(2*n/3)*Delta_x;

first_I = 0 : Delta_x : x0; second_I = x0 + Delta_x : Delta_x : x1; third_I = x1 + Delta_x : Delta_x : L;

IC = zeros(1,n);

a = (x-x1).^2; b = (x-x0).^2; IC = a.*b; IC'

plot(x,IC,'rd')

end

Answer 1

[Cedric]

If IC is a vector, all your expressions

 IC(x) = something

with x a floating point number, are invalid. IC(x) addresses element with index x in vector IC, and valid indices are integers in the range {1,..,length(IC)}. To illustrate:

 >> IC = [7; 8; 9]
 IC =
     7
     8
     9

 >> IC(1)
 ans =
     7

 >> IC(2)
 ans =
     8

 >> IC(3)
 ans =
     9

 >> IC(4)
 Index exceeds matrix dimensions.

 >> IC(0.1)
 Subscript indices must either be real positive integers or logicals.

You can define values for elements of IC based on indies, but I don't think that it was your question. If the question was: you have a vector x and you want to set the initial condition based on values of elements of x, dpb is showing you the way to do it, which I develop in more detail below.

Let's generate a hypothetical x with 10 random elements, and x0 and x1 (I will work with row vectors, because it is easier on the forum):

 >> x = rand(1, 10)
 x =
    0.3404 0.5853 0.2238 0.7513 0.2551 0.5060 0.6991 0.8909 0.9593 0.5472

 >> x0 = 0.3 ;
 >> x1 = 0.7 ;

Now we build a vector of logicals whose true elements will indicate the position of elements in x which are in the range (x0,x1):

 >> id_inRange = x > x0 & x < x1
 id_inRange =
     1     1     0     0     0     1     1     0     0     1

which is the outcome of an AND operation with the two relevant conditions. We will use this vector of logicals to index (logical indexing) relevant elements in both IC and x. Note that you could use pos=find(id_inRange) if you needed to have positions of elements that fall in the range (but logical indexing is better/simpler, so we don't need positions).

Now we want to build IC, first as a vector of zeros (as it is what you seem to want everywhere but where the condition is true):

 >> IC = zeros(size(x))
 IC =
     0     0     0     0     0     0     0     0     0     0

Then, we update relevant elements of IC (those "flagged" by true elements of the logical index), based on corresponding elements of x:

 >> IC(id_inRange) = (x(id_inRange)-x0).^2 .* (x(id_inRange)-x1).^2
 IC =
    0.0002 0.0011  0  0  0  0.0016 0.0000  0  0  0.0014

Hope it helps.

Answer 2

[dpb]

function y=initial_cond(x,x0,x1)
% define function f(x)=(x-x0)^2 * (x-x1)^2 
% over range x0 <= x <= x1;
% f(x) = 0 otherwise.
% INPUT
%    X -- vector of x to be evaluated
%   X0,X1-clipping range
% OUTPUT
%   Y -- vector of length(X) of
%   f(x) = (x-x0)^2 * (x-x1)^2,  x0 <= x <= x1
%        = 0, otherwise

    y=zeros(size(x));
    ix=iswithin(x,x0,x1);
    y(ix)=(x(ix)-x0).^2.*(x(ix)-x1).^2;

This is another use for my utility function iswithin

function flg=iswithin(x,lo,hi)
% returns T for values within range of input
% SYNTAX:
%  [log] = iswithin(x,lo,hi)
%      returns T for x between lo and hi values, inclusive

    flg= (x>=lo) & (x<=hi);

NB: I put the onus of creating the overall x vector and it's spacing to the caller--only the bounds are needed in the function. If you want to pass the start/stop/dx as well, that's a trivial modification obviously.

Answer 3

[dpb]

For the specific inputs so you don't have to think any...

function y = ICondition(L, dx)
  x = 0:dx:L; 
  n = length(x);
  x0 = x(floor(n/3));
  x1 = x(floor(2*n/3));
  ix=iswithin(x,x0,x1);
  y=zeros(size(x));
  y(ix)=(x(ix)-x0).^2.*(x(ix)-x1).^2;

QED

[Eric Johnson]

I appreciate your help, dpb. Thank you.

[dpb]

If you're really going to do it this w/, you might find round better than floor--unless you must have the two intervals <= the value. round will give the closer numerical difference.

[Eric Johnson]

Thanks. I'll try that.

Answer 4

[Eric Johnson]

I can't get that code to run. I need the arguments I had in there because this is a small part of a much larger m-file. Thank you for your help!

[dpb]

Well, mustn't have tried too hard... :(

[Eric Johnson]

The code doesn't run. I get errors... I'm not trying to be rude - it simply fails to run correctly in Matlab.

[dpb]

>> x=[0:0.01:1];
>> x0=1/3;x1=2/3;
>> y=initial_cond(x,x0,x1);
>> plot(x,y)
>> % your original variables instead...
>> 
>> L=1;
>> dx=0.01;
>> x=0:dx:L;
>> y=initial_cond(x,x0,x1);
>> plot(x,y)
>> 
>> % or even closer w/ no changes to the function...
>> y=initial_cond([0:dx:L],x0,x1);
>> plot(x,y)
>> 

function y=initial_cond(x,x0,x1)
% define function f(x)=(x-x0)^2 * (x-x1)^2 
% over range x0 <= x <= x1;
% f(x) = 0 otherwise.
% INPUT
%    X -- vector of x to be evaluated
%   X0,X1-clipping range
% OUTPUT
%   Y -- vector of length(X) of
%   f(x) = (x-x0)^2 * (x-x1)^2,  x0 <= x <= x1
%        = 0, otherwise

    y=zeros(size(x));
    ix=iswithin(x,x0,x1);
    y(ix)=(x(ix)-x0).^2.*(x(ix)-x1).^2;
>>

[dpb]

Sorry if was a little snippy...yesterday had already been a long day at the time so wasn't exactly feeling in a particularly expansive mood... :)

When I had tested that code snippet, hearing that there was some perceived problem in it w/o seeing what had tried to cause the problem I reacted more strongly than perhaps warranted.

Apologies rendered....

Answer 5

[Eric Johnson]

Thank you both for your help. It works!

function[g, IC] = Initial_Condition2(L, Delta_x)

x = 0:Delta_x:L; n = length(x); x0 = floor(n/3)*Delta_x; x1 = floor(2*n/3)*Delta_x; a = (x-x0).^2; b = (x-x1).^2; id_inRange = x > x0 & x < x1 g = a.*b

% you gave me an idea. Multiply the vector g element by element by this id_inRange vector applied to x. Excellent. That should not have taken me 4 hours.

IC = id_inRange.*g

g' IC'

plot(x,g,'b') pause plot(x,IC,'rd')

end

The code I have now isn't exactly beautiful, but it accomplishes what I need it to. (This is a small part of an SSP Runge-Kutta program to test some particle physics problems. The rest I got to work without any issue!) This particular vector is the time derivative for the initial condition U = [1 ... 1]'

[dpb]

Multiply the vector g element by element by this id_inRange vector applied to x.

That's what the logic addressing does transparently...you never showed where you had an error so can only guess but the code posted worked as is. I'd hypothesize the most likely problem is that you didn't use the logical address vector on both LHS and RHS so the effective lengths weren't the same of the target and the object. The above extra multiplication is simply the assignment

g(~idx)=0;

to zero the non-included elements.

The function the way I wrote it just computed the non-zero elements by selecting the subset of the range that is to be included.

[Jan]

@Eric: Please use the "{} Code" button to apply code formatting for a marked section of code. Then your code will be readable.