How to calcul equation using function Matlab

Hi,

I try to write a function for calcul eqution u(x,y t)=exp(-t)*sin(pi*x/L)*sin(pi*y/L) on domaine L-shaped

function z=U(x,y,time,L)
x=[0;0;0;0;0;0;0;0;0;0;0;1;2;3;4;5;6;7;8;9;10;10;10;10;10;10;9;8;7;6;5;5;5;5;5;5;4;3;2;1;1;2;3;4;5;6;7;8;9;1;2;3;4;5;6;7;8;9;1;2;3;4;5;6;7;8;9;1;2;3;4;5;6;7;8;9;1;2;3;4;1;2;3;4;1;2;3;4;1;2;3;4;1;2;3;4];
y=[0;1;2;3;4;5;6;7;8;9;10;0;0;0;0;0;0;0;0;0;0;1;2;3;4;5;5;5;5;5;5;6;7;8;9;10;10;10;10;10;1;1;1;1;1;1;1;1;1;2;2;2;2;2;2;2;2;2;3;3;3;3;3;3;3;3;3;4;4;4;4;4;4;4;4;4;5;5;5;5;6;6;6;6;7;7;7;7;8;8;8;8;9;9;9;9];
z=zeros(size(x));
time=0.15; L=10; E=exp(-time)
%gamma1
 if(x==0 & y>=0&y<=L)
     z=E.*sin(pi*y(x==0&y>=0&y<=L)/L).*sin(pi*0/L);
 
 %gamma2 
elseif( y==L & x>0&x<=L/2)
 z=E.*sin(pi*x(y==L&x>0&x<=L/2)/L).*sin(pi*L/L); 

%gamma3 
elseif( x==L/2 & y>=L/2&y<L)
 z=E.*sin(pi*y(x==L/2&y>=L/2&y<L)/L).*sin(pi*5/L);
 
%gamma4
 elseif( y==L/2 & x>L/2&x<=L)
 z=E.*sin(pi*x(y==L/2&x>L/2&x<=L)/L).*sin(pi*5/L)
 
%gamma5
 elseif (x==L & y>=0&y<L/2) 
z=E.*sin(pi*y(x==L&y>=0&y<L/2)/L).*sin(pi*10/L);
 
 %gamma6 
elseif( y==0 & x>0&x<L)
 z=E.*sin(pi*x(y==0&x>0&x<L)).*sin(pi*0/L);
 
%gamma7
 % elseif ??? I don't how to  write in this case !
z=E.*sin(pi*x/L).*sin(pi*y/L); 
end
end

Please Help me.

Thanks

[dpb]

There is no such array shaped other than as 2D rectangular array so you have no way to store your z variable as a single array in such a form. You would have to pick a storage order and index in some fashion or fill in the missing corner values with NaN or somesuch.

The simplest way to do what you're asking would be to use meshgrid and evaluate the function over the entire outer boundaries first, then reset the upper quadrant to NaN to indicate the missing locations.

If you don't want to return those in that fashion, then at that point you could make the storage transformation to whatever holding arrangement you did choose.

What is the intended use of the result? That could help in suggesting storage alternatives.

[Yamina chbak]

Yes, you are right. I use delaunary triangulation and I have Nodes, elements and edges for use triplot to make a simple tringular mesh of the L-shaped. Here i give you a list of x,y and TRI to see resultat

x=[0;0;0;0;0;0;0;0;0;0;0;1;2;3;4;5;6;7;8;9;10;10;10;10;10;10;9;8;7;6;5;5;5;5;5;5;4;3;2;1;1;2;3;4;5;6;7;8;9;1;2;3;4;5;6;7;8;9;1;2;3;4;5;6;7;8;9;1;2;3;4;5;6;7;8;9;1;2;3;4;1;2;3;4;1;2;3;4;1;2;3;4;1;2;3;4];
y=[0;1;2;3;4;5;6;7;8;9;10;0;0;0;0;0;0;0;0;0;0;1;2;3;4;5;5;5;5;5;5;6;7;8;9;10;10;10;10;10;1;1;1;1;1;1;1;1;1;2;2;2;2;2;2;2;2;2;3;3;3;3;3;3;3;3;3;4;4;4;4;4;4;4;4;4;5;5;5;5;6;6;6;6;7;7;7;7;8;8;8;8;9;9;9;9];
TRI=[1 12 2;12 2 41;12 13 42;12 42 41;13 14 42;14 43 42;14 43 44;14 15 44;15 44 16;44 16 45;45 46 16;16 17 46;46 18 17;46 47 18;18 48 47;18 19 48;19 20 48;48 20 49;20 22 49;20 21 22;2 3 50;2 41 50;41 42 50;42 50 51;42 51 52;42 43 52;44 52 53;43 44 52;44 54 53;44 45 54;45 46 54;46 55 54;46 47 56;46 55 56;47 48 56;48 57 56;57 58 48;48 49 58;49 22 58;23 22 58;50 4 59;3 4 50;59 60 50;50 51 60;61 52 60;51 52 60;53 52 62;52 61 62;62 63 54;53 54 62;54 64 63;54 55 64;56 65 64;55 56 64;56 66 65;56 57 66;66 58 67;57 58 66;58 24 67;58 23 24;4 68 5;4 59 68;68 69 60;59 68 60;60 70 69;60 61 70;62 71 70;61 70 62;72 62 71;62 63 72;64 73 72;63 64 72;64 74 73;64 65 74;66 75 74;65 66 74;66 76 75;66 67 76;24 76 25;24 76 67;68 6 5;68 6 77;77 78 68;69 78 68;69 70 78;79 70 78;71 80 70;79 80 70;71 72 80;31 72 80;72 30 73;72 30 31;74 30 73;74 30 29;74 28 75;74 28 29;76 28 75;76 28 27;76 26 25;76 26 27;6 81 77;6 81 7;78 81 77;78 81 82;78 83 79;78 83 82;80 83 79;80 83 84;80 32 31;80 32 84;81 8 7;81 8 85;81 86 82;81 86 85;83 86 82;83 86 87;83 88 84;83 88 87;32 88 84;32 88 33;89 8 85;89 8 9;85 89 86;90 89 86;86 91 87;86 91 90;91 88 87;91 88 91;88 34 33;88 34 92;9 10 89;93 10 89;94 89 90;94 89 93;90 94 91;95 94 91;92 91 96;95 91 96;34 96 92;34 96 35;40 10 93;40 10 11;94 40 93;94 40 39;95 38 94;39 38 94;38 96 95;38 96 37;96 36 35;96 36 37];
figure(1)
triplot (TRI, x, y)

But, I sill i don't khow how to calcul equation U on L-shaped domain. (It's about boundaries condition and initial conditions)

[Yamina chbak]

@dpb, I have modified it and please check it out and help me. Thanks.

[John D'Errico]

There should be NO problem. Just compute it on the entire domain, and then set the elements outside the L-shaped boundary to NaN. So where is the disconnect?

Can you compute the result on the entire domain? A simple vectorized computation would do that, in as little as one line of code. Then, can you set the elements to NaN in a corner?

[Yamina chbak]

@John D'errico, Did you mean compute it on, for example, square domain ? And set the elements outside the L-shaped to NaN? But this takes me written a two mesh of the square and the L-shaped.

In fact, My goal is to compute U on the L-shaped domain with a large number of elements for use it my other code.

The most important is to write a two function: first, compute U on the L-shaped boundary. When I call it, I find Error: Output argument "z" not assigned durig call ....

Second, If it is working ( I hope it), I make another function to compute it on the L-shaped (with gamma7, see my figure above). Here I don't know how to write in this case (gamma7)!

Please help me.

[dpb]

"... you mean compute it on, for example, square domain ? And set the elements outside the L-shaped to NaN?"

That's precisely what we meant, yes. There's certainly no easier way to compute the complete domain and you can extract the boundary segment values from it if returning those is another desire.

Those alone could be computed far more efficiently by just computing along the specific lines without all the logical addressing for whatever spacing along the x,y axes is desired.

Your function as written doesn't have any way to return but the last z vector because each clause in the if...elseif...end block writes to the same variable.

The if block isn't properly constructed; in MATLAB an IF is satisfied IFF all elements of the expression are true; your x and y being arrays will return an array of logicals only a few of which will be T. Hence, none of the expressions will be satisfied.

Lastly, there's no point in having x,y,time,L as arguments in the function since you immediately redefine them inside the function.

[Yamina chbak]

I get your point, I will try it if i have no choice, I can compute the resultat on the square domain but i don't know how to set the elements outside the L-shaped to NaN in a corner!. I am a beginner at Matlab.

Answer 1

[Walter Roberson]

function z=U(x,y,time,L)
x=[0;0;0;0;0;0;0;0;0;0;0;1;2;3;4;5;6;7;8;9;10;10;10;10;10;10;9;8;7;6;5;5;5;5;5;5;4;3;2;1;1;2;3;4;5;6;7;8;9;1;2;3;4;5;6;7;8;9;1;2;3;4;5;6;7;8;9;1;2;3;4;5;6;7;8;9;1;2;3;4;1;2;3;4;1;2;3;4;1;2;3;4;1;2;3;4];
y=[0;1;2;3;4;5;6;7;8;9;10;0;0;0;0;0;0;0;0;0;0;1;2;3;4;5;5;5;5;5;5;6;7;8;9;10;10;10;10;10;1;1;1;1;1;1;1;1;1;2;2;2;2;2;2;2;2;2;3;3;3;3;3;3;3;3;3;4;4;4;4;4;4;4;4;4;5;5;5;5;6;6;6;6;7;7;7;7;8;8;8;8;9;9;9;9];

z=zeros(size(x));

time=0.15; L=10; E=exp(-time)

%gamma1
 mask = (x==0 & y>=0&y<=L)
 z(mask) = E.*sin(pi*y(mask)/L).*sin(pi*0/L);
 
 %gamma2 
mask = ( y==L & x>0&x<=L/2)
 z(mask) = E.*sin(pi*x(mask)/L).*sin(pi*L/L); 

%gamma3 
mask = ( x==L/2 & y>=L/2&y<L)
z = E.*sin(pi*y(mask)/L).*sin(pi*5/L);

and so on.

[Yamina chbak]

@Walter Roberson, What about the case of gamma7 (see my figure above), I mean is computed the equation inside an L-shaped domain, to get exactly that size(z)=size(x).

[Walter Roberson]

function z=U(x,y,time,L)
x=[0;0;0;0;0;0;0;0;0;0;0;1;2;3;4;5;6;7;8;9;10;10;10;10;10;10;9;8;7;6;5;5;5;5;5;5;4;3;2;1;1;2;3;4;5;6;7;8;9;1;2;3;4;5;6;7;8;9;1;2;3;4;5;6;7;8;9;1;2;3;4;5;6;7;8;9;1;2;3;4;1;2;3;4;1;2;3;4;1;2;3;4;1;2;3;4];
y=[0;1;2;3;4;5;6;7;8;9;10;0;0;0;0;0;0;0;0;0;0;1;2;3;4;5;5;5;5;5;5;6;7;8;9;10;10;10;10;10;1;1;1;1;1;1;1;1;1;2;2;2;2;2;2;2;2;2;3;3;3;3;3;3;3;3;3;4;4;4;4;4;4;4;4;4;5;5;5;5;6;6;6;6;7;7;7;7;8;8;8;8;9;9;9;9];

z=zeros(size(x));

time=0.15; L=10; E=exp(-time)

%gamma1
 mask1 = (x==0 & y>=0&y<=L)
 z(mask1) = E.*sin(pi*y(mask1)/L).*sin(pi*0/L);
 
 %gamma2 
mask2 = ( y==L & x>0&x<=L/2)
 z(mask2) = E.*sin(pi*x(mask2)/L).*sin(pi*L/L); 

%gamma3 
mask3 = ( x==L/2 & y>=L/2&y<L)
z = E.*sin(pi*y(mask3)/L).*sin(pi*5/L);

etc

mask7 = ~(mask1 | mask2 | mask3 | mask4 | mask5 | mask6);
z(mask7) = E.*sin(pi*x(mask7)/L).*sin(pi*y(mask7)/L); 

[Yamina chbak]

Thanks you @Walter Roberson. It's working.