## Posts Tagged ‘Matlab

21
Oct
06

### Vortices

Here is another code to solve Euler’s equation (inviscid flow) and plot the pathlines. This method was adapted from Stephen Montgomery-Smith’sEuler2D program. The vortex effects are found using the Biot-Savart law and the differential equations are solved using the Adams Bashforth method. Click here to download some Matlab code for this. Shown below is some Mathematica code:
```(* runtime: 19 seconds *) ToPoint[z_] := {Re[z], Im[z]};dt = 0.01; rcore = 0.1; dz0 = 0; Klist = {1, -1, 1, -1}; zlist = Join[{-1 - 0.5I, -1 + 0.5I, -0.5 - 0.5I, -0.5 + 0.5I}, Flatten[Table[x + I y, {x, -1, 1, 0.1}, {y, -1, 1, 0.1}], 1]]; paths = Map[# Table[1, {10}] &, zlist]; Do[dz = Map[Sum[zj = zlist[[j]]; If[# != zj, r2 = Abs[# - zj]^2; I Klist[[j]](zj - #)/r2 (1 -Exp[-r2/rcore^2]), 0], {j, 1, Length[Klist]}] &, zlist]; zlist += dt (1.5dz - 0.5dz0); dz0 = dz; paths = Table[Prepend[Delete[paths[[i]], -1], zlist[[i]]], {i, 1, Length[zlist]}];Show[Graphics[Map[Line[Map[ToPoint, #]] &, paths], PlotRange -> {{-1, 1}, {-1, 1}}, AspectRatio -> Automatic]], {85}];```

12
Jun
06

### Unstructured Tetrahedral Meshed Sphere

Here is an unstructured tetrahedral meshed sphere.

24
May
06

### Finite Element Method (FEM) Solution to Poisson’s equation on Triangular Mesh

Finite Element Method (FEM) Solution to Poisson’s equation on Triangular Mesh
solved in Mathematica 4.2, 5/24/06; mesh generated with Gmsh; old version: Matlab and DistMesh

Here is some code to solve Poisson’s equation on an unstructured grid of triangular elements using the Finite Element Method (FEM):
```(* runtime: 0.1 second *) n = 10; nodes = Flatten[Table[{x, y}, {y, -1, 1, 2.0/(n - 1)}, {x, -1, 1, 2.0/(n - 1)}], 1]; elements = Flatten[Table[{{i, i + 1, i + n}, {i + 1, i + n + 1, i + n}} + (j - 1) n, {j, 1, n - 1}, {i, 1, n - 1}], 2]; nnode = Length[nodes]; nelem = Length[elements]; fixed = Map[Position[nodes, #][[1, 1]] &, Select[nodes, Max[Abs[#]] == 1 &]]; Kglobal = Table[0, {nnode}, {nnode}]; F = phi = Table[0, {nnode}]; Scan[(xy =Join[Transpose[nodes[[#]]], {{1, 1, 1}}];deta = Inverse[xy][[All, {1, 2}]]; Kglobal[[#, #]] +=Det[xy]deta.Transpose[deta]/2; F[[#]] += Det[xy]/6) &, elements]; free = Complement[Range[nnode], fixed]; phi[[free]] = LinearSolve[Kglobal[[free, free]], F[[free]]]; Mean[x_] := Plus @@ x/Length[x]; PlotColor[x_] := Hue[2(1 - x)/3]; Show[Graphics[Table[{PlotColor[Mean[phi[[elements[[i]]]]]/Max[phi]], Polygon[nodes[[elements[[i]]]]]}, {i,1, nelem}], AspectRatio -> 1]];```

Here is some code for an interpolated plot:
```(* runtime: 10 seconds *) n = 275; x1 = y1 = -1; x2 = y2 = 1; image = Table[0, {n}, {n}]; xIntersect[{{x1_, y1_}, {x2_, y2_}}] := If[y1 == y2, Infinity, x1 + (y - y1)(x2 - x1)/(y2 - y1)]; Scan[(plist = nodes[[#]]; xlist = nodes[[#, 1]]; ylist = nodes[[#, 2]]; p1 = plist[[1]]; {dx2, dy2} = plist[[2]] - p1; {dx3, dy3} = plist[[3]] - p1; Do[y = y1 + (y2 - y1)(i - 1)/(n - 1); {xa, xb,xc} = Map[xIntersect, Partition[Sort[plist, #1[[2]] < #2[[2]] &], 2, 1, 1]]; jlist = Round[n({xc, If[Abs[xb - xc] < Abs[xa - xc], xb, xa]} - x1)/(x2 - x1)]; If[jlist =!= {}, Do[x = x1 + (x2 - x1)(j - 1)/(n - 1); {xi, eta} = ((x - p1[[1]]){dy3, -dy2} + (y - p1[[2]]){-dx3, dx2})/(dx2 dy3 - dx3 dy2); image[[i, j]] = phi[[#]].{1 - xi - eta, xi, eta}, {j, Max[1, Min[jlist]], Min[n, Max[jlist]]}]], {i, Max[1, Round[n(Min[ylist] - y1)/(y2 - y1)]], Min[n, Round[n(Max[ylist] - y1)/(y2 - y1)]]}]) &, elements]; ListDensityPlot[image, Mesh -> False, Frame -> False, ColorFunction -> PlotColor, Epilog -> Map[{Line[Map[({x, y} = nodes[[#, {1, 2}]]; n{(x - x1)/(x2 -x1), (y - y1)/(y2 - y1)}) &, #]]} &, elements]];```

## Welcome !

You will find here some of my favorite hobbies and interests, especially science and art.

I hope you enjoy it!

Subscribe to the RSS feed to stay informed when I publish something new here.

I would love to hear from you! Please feel free to send me an email : bugman123-at-gmail-dot-com