直径3mmの竹ひごの曲げを調べる 竹ひごの長さは11.8mm Matlabを使ってピクセルでの距離を求め、自然長からたわみを求める。 grams = [0 60 80]; ys = []; for i=1:length(grams) fname = string(grams(i))+'gf.png'; im = imread(fname); imshow(im); [x,y,b…

スキンをサンドイッチ構造で作る www.jcrocket.com

ストリンガーで補強したのがセミモノコック https://www.google.com/url?sa=t&rct=j&q=&esrc=s&source=web&cd=&ved=2ahUKEwilssK9ysvrAhWTZt4KHREQBkUQFjAAegQIAxAB&url=http%3A%2F%2Fwww.aspirespace.org.uk%2Fdownloads%2FRocket%2520airframe%2520concept…

LU分解 →Aが正則行列のとき 修正コレスキー分解 →Aがエルミート行列のとき コレスキー分解 →Aがエルミート行列かつ正定値行列のとき ここで、行列[M_k]は以下で定義する。 →ピボットで割る計算例 clear; A = [10 6 2; 6 5 2; 2 2 1]; b = [1 2 3]'; M1 = [1 …

close all; clear; n=500; r1s=[]; r2s=[]; xs=[]; ts=1:n; for t=1:n r1=rand(); r2=rand(); x = sqrt(-2*log(r1))*sin(2*pi*r2); r1s(end+1)=r1; r2s(end+1)=r2; xs(end+1)=x; end subplot(211) plot(ts,r1s); hold on; plot(ts,r2s); plot(ts,xs); legend…

非線形自由振動 初期条件 とする解析解 clear; close; t=0:0.1:10 omega=1 e=0.1 for a=[1,5,10,11,12] x = a*cos(omega*t)+... e*a^3/(32*omega^2)*(cos(3*omega*t)-cos(omega*t))+... e^2*a^5./(1024*omega^4)*(cos(5*omega*t)-cos(omega*t))+... e^3*a^7/…

studentcommunity.ansys.comANSYS WorkbenchのOptionでBeta optionsを有効にする SymmetryからGraphical expansion(Beta)のNum Repeatを２に、MethodをHalfに、⊿Zを1e-3に設定してみたが、表示結果は変わらなかった。

clear; close all; figure();grid(); xlim([-10 10]);ylim([-10 10]); [x,y]=ginput(4); x0=x(1);x1=x(2);x2=x(3);x3=x(4); y0=y(1);y1=y(2);y2=y(3);y3=y(4); Px1=[];Py1=[]; Px2=[];Py2=[]; for t=0:0.01:1 %Bézier B0=(1-t)^3; B1=3*(1-t)^2*t; B2=3*(1-t…

clear all; close all; q01=2; q02=-2; e0=8.85*1e-12; k=1/(4*pi*e0); omega=1; r1=[-2,0]'; r2=[2,0]'; x=-5:5; y=-5:5; [X,Y]=meshgrid(x,y); a1=sqrt((X-r1(1)).^2+(Y-r1(2)).^2); a2=sqrt((X-r2(1)).^2+(Y-r2(2)).^2); myVideo = VideoWriter('myVideoF…

clear all; close all; q1=2; q2=-2; e0=8.85*1e-12; k=1/(4*pi*e0); x=-5:5; y=-5:5; r1=[-2,0]'; r2=[2,0]'; [X,Y]=meshgrid(x,y); a1=sqrt((X-r1(1)).^2+(Y-r1(2)).^2); a2=sqrt((X-r2(1)).^2+(Y-r2(2)).^2); E1x=q1*k./a1^3.*(X-r1(1)); E1y=q1*k./a1^3.…

clear all; close all; I=1; x=-5:5; y=-5:5; z=-5:5; [X,Y,Z]=meshgrid(x,y,z); Hx=-I/(2*pi)*Y./(X.^2+Y.^2); Hy=I/(2*pi)*X./(X.^2+Y.^2); Hz=0*X; plot3([0,0],[0,0],[z(1),z(end)],'r-','linewidth',5); hold on; quiver3(X,Y,Z,Hx,Hy,Hz,'b');

clear all; close all; % q1=2; % q2=-2; e0=8.85*1e-12; k=1/(4*pi*e0); q1s=[1,2,3]; q2s=[-3,-2,-1]; x=-5:5; y=-5:5; r1=[-2,0]'; r2=[2,0]' [X,Y]=meshgrid(x,y); a1=sqrt((X-r1(1)).^2+(Y-r1(2)).^2); a2=sqrt((X-r2(1)).^2+(Y-r2(2)).^2); for i1=1:l…

clear all; close all; q=2; x=-10:10; y=-10:10; [X,Y]=meshgrid(x,y); Z=sqrt(100-X.^2-Y.^2); Ex=100*X./(X.^2+Y.^2+Z.^2).^(3/2); Ey=100*Y./(X.^2+Y.^2+Z.^2).^(3/2); Ez=100*Z./(X.^2+Y.^2+Z.^2).^(3/2); Ex(isnan(Ex))=0; Ey(isnan(Ey))=0; Ez(isnan(…

clear all; close all; q=2; r=[[1,0]',[0,1]',[-1,0]',[0,-1]',[1.5,1.5]',[-1.5,1.5]',[-1.5,-1.5]',[1.5,-1.5]']; r2=arrayfun(@(x,y)x.^2+y.^2,r(1,:),r(2,:)) r1=sqrt(r2); e0=8.85*1e-12; E=1/(4*pi*e0)*q./r2./r1.*r; plot(0,0,'bo'); hold on; array…

clear all; close all; qA=-2; qB=1; rA=[2,1,0]'; rB=[0,1,0]'; r=[rA rB]; e0=8.85*1e-12; k = 1/4/pi/e0; F = k*qA*qB/sqrt(sum((rA-rB).^2)); Fdummy=1;%correction for view F=Fdummy;%The real F is too huge to be shown, so just leave it v=rA-rB; …

よさそう板の曲げと静荷重弾性解析ウィザード https://sites.google.com/site/codeastersalomemeca/home/salome-meca-chutoriaru1 eficas（コマンドファイルエディタ）の使用方法 https://sites.google.com/site/codeastersalomemeca/home/salome-meca-chuto…

ヤコビ法 function jacobi % dim = 10; % A = reshape(rand(dim),[dim,dim]); % b = rand(dim,1); A= [10 3 1 2 1 1 19 2 -1 5 -1 1 30 1 10 -2 0 1 20 5 -3 5 1 -2 25]; b = [-22 27 89 -73 22]'; % x = A\b [L,D,U] = decomposeLDU(A); % testLDU(L,D,U); …

Conjugate Gradient(CG: 共役勾配法) Hierarchical Domain Decomposition Method (HDDM: 階層型領域分割法) Balancing Domain Decomposition (BDD前処理)

https://github.com/lukeasrodgers/frepl

なかなかすごい需要だ https://github.com/jacobwilliams/pyplot-fortranFoBiSっていうFortranのビルドシステムがあることを知ったhttps://github.com/szaghi/FoBiS program test use,intrinsic :: iso_fortran_env, only: wp => real64 use pyplot_module i…

http://nineties.github.io/math-seminar/6.html#/31

・内挿＝補間＝interpolation ・外挿＝補外＝extrapolation

function chebyshev figure();hold on; title('Chebyshev polynomial'); for n=1:4 theta = linspace(0,2*pi,100); x = cos(theta); Tnx = cos(n*theta); plot(x,Tnx); end l = legend({'N=1' 'N=2' 'N=3' 'N=4'}); set(l,'FontSize',14); end チェビシェフ…

結構めんどくさい うまく動いてない途中で三重対角行列を係数行列とする連立１次方程式の解を求める必要があったが、TDMA（三重対角行列アルゴリズム）という解法がある。 function spline_interpolation clear all; close all; format compact; div = 100; …

N=3ですでにサインカーブがうまく補間されている function lagrange_interpolation clear all; close all; div = 100; x = linspace(0,pi,100); y = sin(x); plot(x,y); hold on; a = x(1); b = x(end); fx = (x-b)/(a-b)*sin(a) + (x-a)/(b-a)*sin(b); plot…

http://www.takuichi.net/hobby/edu/em/mom/numerical/index-j.htm を計算する リーマン積分（区分求積法） function riemann clear all; close all; divlist = [10e2,10e3,10e4,10e5,10e6,10e7]; errorlist = []; for cycle=1:6 div = divlist(cycle); x = …

function conway2d clear all; close all; CAsize = 100; [X,Y] = meshgrid(1,CAsize+1); figure(); hold on; plot(X,Y,'k'); plot(Y,X,'k'); I = zeros(CAsize+1); surface(I); colormap(gray); axis off; I = rand(CAsize+1)>0.85; I(1,:) = 0; I(CAsize,:…

1-256の規則を全て適用してみた function wolfram1d_2 clear all; close all; CAsize = 100; [X,Y] = meshgrid(1,CAsize+1); k = {'111','110','101','100','011','010','001','000'}; % v = {'1','0','1','0','0','1','0','1'}; %165 % writerObj = VideoWr…

規則182を適用 function wolfram1d_2 clear all; close all; CAsize = 100; [X,Y] = meshgrid(1,CAsize+1); figure(); hold on; plot(X,Y,'k'); plot(Y,X,'k'); I = zeros(CAsize+1); surface(I); colormap(gray); axis off; origin = zeros(1,CAsize); % fo…

function velocity_variable_diffusion clear all; close all; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% CAPTURE = false; VIEWTOP = true; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %condition …