Matlab Thread

[Anonymous No.5112969]

Post and discuss cool Matlab scripts. Include short description at top of code posts. If code is too long use pastebin? OC only.

I'll start slow, here is a 1D wave on a string you can play with your mouse:

clear all, close all; clc

c = 5; % wave speed
dx = .01;
dt = .01;
n = 101;
len = dx*(n - 1);
x = 0:dx:len; % x dim for plotting
y = 768; % screen pixel height
r2 = (c*dt/dx)^2;

Ko = diag(ones(1,n-1),-1) + diag(-2*ones(1,n),0) + diag(ones(1,n-1),1);
K = sparse(Ko*r2);
K([1,end],:) = 0;
A = (speye(n,n) - K);

uold = zeros(n,1);
ucurr = uold;
t = 0;

figure('Units','normalized','Position',[.1 .1 .8 .8])
while t <= 100
mxy = get(0,'PointerLocation');
ucurr(1) = 2*mxy(2)/y - 1;
unew = A\(2*ucurr - uold);
uold = ucurr;
ucurr = unew;
t = t + dt;

plot(x(1:end-1),unew(2:end))
axis([0 len -1 1]), getframe;
end

[Anonymous No.5113036]

>>5112969
That was awesome.

[Anonymous No.5113056]

That was pretty cool OP.

[Anonymous No.5113697]

here's another one, "an implicit meshing algorithm".Please excuse the lack of vectorization. Also, I am aware that all the forces are being computed twice, I have versions that run faster and on any 2D boundary, this is just the compact version. Enjoy the animation!

clear all, close all; clc
figure('units','normalized','position',[.2 .2 .6 .6])

r = .5; % Boundary radius
k = -.5; % Spring constant
h = .5; % 'Time' step
len = .08; % Spring element rest length
N = 200; % Number of nodes
A = rand(N,2); % Generating N random x-y points [0 to 1]
A = (A*(r*2) - r)/2; % Scaling and shifting points
for it = 1:75
fx = zeros(1,N);
fy = zeros(1,N);
Tri = delaunay(A(:,1),(A(:,2)));
for n = 1:N
for i = 1:length(Tri)
for j = 1:3
if Tri(i,j)==n
l = [1,2,3,1,2];
lenx1 = (A(Tri(i,l(j+1)),1) - A(Tri(i,j),1));
lenx2 = (A(Tri(i,l(j+2)),1) - A(Tri(i,j),1));
leny1 = (A(Tri(i,l(j+1)),2) - A(Tri(i,j),2));
leny2 = (A(Tri(i,l(j+2)),2) - A(Tri(i,j),2));
nrm1 = (lenx1^2 + leny1^2)^.5;
nrm2 = (lenx2^2 + leny2^2)^.5;
fx(n) = fx(n)+k*((len-nrm1)*(lenx1/nrm1) + (len-nrm2)*(lenx2/nrm2));
fy(n) = fy(n)+k*((len-nrm1)*(leny1/nrm1) + (len-nrm2)*(leny2/nrm2));
end
end
end
A(n,1)= A(n,1)+h*fx(n);
A(n,2)= A(n,2)+h*fy(n);
end
for n = 1:N
NRM = (A(n,1)^2 + A(n,2)^2)^.5;
if NRM > r
A(n,:) = A(n,:)*r/NRM;
end
end
triplot(Tri,A(:,1),(A(:,2))), hold on
[x,y,z] = cylinder(r,200);
plot(x(1,:),y(1,:),'r')
axis equal, hold off
title('Implicit Meshing of a Circle')
getframe;
end

[Anonymous No.5113749]

>>5112969
that's pretty cool OP. like to see more.

[Anonymous No.5113784]

someone have more? :)

[Anonymous No.5113794]

>>5112969
>Using \ correctly

you get an A in my book

- A Concerned MathWorks employee

[Anonymous No.5113799]

>>5113794
how did i use \ incorrectly?

[Anonymous No.5113802]

error: get: unknown property "PointerLocation"

;<

[Anonymous No.5113808]

>>5113794
Wait, so what does it do exactly?

[Anonymous No.5113810]

>>5113799
haha nvm, i misread, here is 2D version of original post, you will like this one even moar!

clear all, close all
% ================= problem params: =====================
c = 7; % c = [m/s]
dx = .1; % dx = dy
m = 200; % x nodes
n = 200; % y nodes
Lx = (m-1)*dx; % x length
Ly = (n-1)*dx; % y length
x = 0:dx:Lx;
y = 0:dx:Ly;
to = 0;
tf = 1000;
dt = .01;
t = 0:dt:tf;
r2 = (dt*c/dx)^2;
v = r2*[1,1,-4,1,1]; % vector of coeffs

d = [-n -1 0 1 n];
B = ones(m*n,1)*v;
Ks = spdiags(B,d,m*n,m*n);
ij2node = reshape((1:m*n),n,m)';
lbn = ij2node(:,1);
rbn = ij2node(:,end);
tbn = ij2node(1,:)';
bbn = ij2node(end,:)';
bn = unique([lbn;rbn;tbn;bbn]);

K = (Ks + 2*speye(m*n,m*n));

u = 0*ones(m*n,1);
H = .2;
Uold = u;
Ucurr = Uold;
figure('Units','normalized','Position',[.25 .25 .5 .5])
% =================== solve: ===========================
for tt = 2:length(t)-1

mxy = get(0,'PointerLocation');
ni = ij2node(floor(3 + mxy(1)*(m)/1920),(3 + floor(mxy(2)*(n)/1080)));
Ucurr(ni+d) = H*.9;
Ucurr(ni) = H;
Ucurr(bn) = 0;
Unew = K*Ucurr - Uold;
Unew(ni+d) = H*.9;
Unew(ni) = H;
Uold = Ucurr;
Ucurr = Unew;
Unew(end) = 1;
Unew(1) = -.3;

sol = reshape(Unew,n,m);
surf(x,y,sol), shading interp
axis equal, zlim([-2 2])
getframe;
end

[Anonymous No.5113820]

>>5113808

backslash is a linear system solving operator, it can use a variety of methods depending on the structure of the system (banded, dense, sparse, etc...)

[Anonymous No.5113827]

>>5113799
no no, there was no sarcasm in my post. most people don't know how to use it and i am glad to see that someone does.

[Anonymous No.5113837]

>>5113810
holy shit! awesome! i would've never thought of doing stuff like this in matlab.

[Anonymous No.5113843]

>>5113802
you must have some old version of matlab or something. sorry bro, try the second code instead, no pointer location needed.

[Anonymous No.5113862]

>>5113820
WELL FUCK it looks useful.
Why don't you make it more obvious, I never stumbled upon it in the help.
Obviously, I'm not going to do a search for "\" if I don't already know it exists.

[Anonymous No.5113866]

>>5113820
Wait you can use proper banded structures in Matlab like in Lapack?
How come I never found it and had to do a fucking Fortran script then?

[Anonymous No.5113883]

>>5113866
check out this mathworks page:

http://www.mathworks.com/help/matlab/ref/spparms.html

[Anonymous No.5114011]

>>5113883
Oh, that kind of band.
I was thinking of dense band matrices like in ZHBEV

[Anonymous No.5114271]

OP here again. I really want someone else to post some code!

Here is a cute little n-body gravity simulation that I know you will love:

clear all, close all; clc
figure('units','normalized','position',[.1 .1 .7 .7])

to = 0;
tf = 10;
dt = .01;
T = to:dt:tf;
n = 100;
a = 1e3*.2*dt;
r = 20;

Xold1 = [r*2*(rand(n,3) - .5)] + 50;
Xold2 = [r*2*(rand(n,3) - .5)] - 50;

Xold = [Xold1;Xold2;0 0 0];
Xc = Xold*[cos(a*pi/180) -sin(a*pi/180) 0;
sin(a*pi/180) cos(a*pi/180) 0;
0 0 1];

M = [1e2*(rand(2*n,1).^5); 1e7];

N = length(M);
dt2 = dt^2;
Xn = zeros(N,3);

for ti = 1:length(T)
F = zeros(N,3);
for n = 1:N
I = [1:n-1,n+1:N];
for i = 1:N-1
V = Xc(I(i),:) - Xc(n,:);
norm_V = norm(V);
F(n,:) = F(n,:) + (M(I(i))/(norm_V^3))*V;
end
end
Xn = 2*Xc - Xold + dt2*F;
Xold = Xc;
Xc = Xn;
plot3(Xn(:,1),Xn(:,2),Xn(:,3),'k*'), grid on
axis([-150 150 -150 150 -150 150]), getframe;
end

[Anonymous No.5114294]

>>5114271
Aw dude that's so cool.

[Anonymous No.5114380]

I have a ~4000 line MRI segmentation algorithm

>implying I'm going to post that fucker

[Anonymous No.5114393]

>>5114380

post something cute with cool visuals, you dont see me posting my fast multipole accelerated boundary element code...

[Anonymous No.5115522]

>>5114271
bump
I love your elegant coding OP.

[Anonymous No.5115547]

babby's first MATLAB here, was trying to make a least squares function

what's the difference between

beta = (x'*x)^-1 *(x' * y);
and
beta = (x'*x) \ (x' * y);

?

[Anonymous No.5115586]

>>5115547
\ is pretty much always more efficient because actually calculating the inverse is really slow. Use backslash.

[Anonymous No.5115755]

Give me $1000 for the licensing and then I'll figure out some coold codes OP.

[Anonymous No.5115775]

Anyone got any MATLAB tutorials?
Incredibly new to MATLAB and need to learn how to use it.
Like, how to do the Netwon-Raphson method with it.

[Anonymous No.5115793]

>>5115547
uhm, (x'*x)^-1 *(x' * y) <=> (x' * y) / (x' * x) =/= (x'*x) \ (x' * y)

[Anonymous No.5115838]

>actually using proprietary software

Non-free software is an enemy of intellectual freedom.

I recommend GNU Octave, a high level language similar to matlab that respects your essential freedoms.

>https://www.gnu.org/software/octave/

[Anonymous No.5115840]

>>5115838
>GNU Octave
>not Scilab

[Anonymous No.5115848]

>>5115840

Or this, thank you based anon

[Anonymous No.5115908]

>>5115840
>>5115848

> Be same t-shirt
> Picture taken ten years apart
> Never washed.

[Anonymous No.5115917]

>>5115908
In my land, we call that a polo.

[Anonymous No.5115944]

Hey OP, love the thread premise and all, but have to comment
>my exorbitantly disfigured face when I'm uncontrollably laughing at people who have a tough time in Matlab intro courses
I realize I've done Java I already and have a lot of Lua experience, but still it isn't that hard if you give the class an iota of your efforts. The teacher is pretty terrible though, a bad way to learn programming.

[Anonymous No.5115950]

>>5115917
It's not a polo.
It's just torn and has mold that makes it appear as such

[Anonymous No.5117964]

last call for anyone else to contribute

[Anonymous No.5118015]

>>5117964
I WISH I COULD

[Anonymous No.5118073]

>>5114271

Frickin awesome! Is there a supermassive center?

I just installed matlab today. Im saving your code so I can learn from it :)

Please enjoy a glance at this log plot of the 26 heaviest known bodies in our solar system.

[Anonymous No.5118141]

I used MATLAB for some intro math courses, and haven't touched it since. have used Python quite a lot though and I like it. I'm not doing anything so intense that the slow run times are a problem.

[Anonymous No.5118161]

what was the reasoning behind matlab showing every result from a line of code unless you inlcude a ' ; ' on the end? Who the fuck thought that was a good idea?

[Anonymous No.5118354]

>>5118073

yeah

> M = [1e2*(rand(2*n,1).^5); 1e7];

1e7 is the massive center and the other term is "n" random masses. raising to the ".5" shifts the distribution if you think about it.

[Anonymous No.5118417]

>>5114271
pretty cool

[Anonymous No.5118447]

I can't find the code now, but one of the assignments a few years back I had was to make an analogue clock that told the current time using MATLAB (by using plots of circles and lines etc).

[Anonymous No.5118660]

Anyone got any good guides on how to use MATLAB?