以下のローレンツモデルを解く clear; close all; a=10; b=8/3; r=28; f1 = @(t,x,y,z) -a*x+a*y; f2 = @(t,x,y,z) -x*z+r*x-y; f3 = @(t,x,y,z) x*y - b*z; N = 2000; t0 = 0; tN = 40; h = (tN-t0)/N; t = linspace(t0,tN,N); x = zeros(1,N); y = zeros(1…

を計算する。 解析解は,rungekutta.m clear; close all; x0=0; xN=1; N=10; h=(xN-x0)/N; x=linspace(x0,xN,N); y=zeros(1,N); a = 1; y(1)=1; f = @(x,y) a*y; for k=2:N k1 = h*f(x(k-1),y(k-1)); k2 = h*f(x(k-1)+h/2,y(k-1)+k1/2); k3 = h*f(x(k-1)+h/2,…

を解く。 解析解は、 clear; close all; x0=0; xN=1; N=20; h=(xN-x0)/N; x=linspace(x0,xN,N); y=zeros(1,N); %syms y0; a = 1; y(1)=1; f = @(x,y) a*y; for k=2:N k1 = h*f(x(k-1),y(k-1)); k2 = h*f(x(k-1)+h,y(k-1)+k1); y(k) = y(k-1) + 1/2*(k1+k2); …

を解く。 初期条件はおよび 上記は以下の連立微分方程式に変形できる。 clear all; close all; %dydx = @(x,y) -y -pi*exp(-x)*sin(pi*x); N = 20; x0 = 0; xN = 10; %%%% Forward Euler %%%% x = linspace(x0,xN,2*N); y = zeros(1,2*N); z = zeros(1,2*N);…

を解く。 この方程式の解析解は、 clear all; close all; dydx = @(x,y) -y -pi*exp(-x)*sin(pi*x); N = 20; x0 = 0; xN = 2; %%%% Forward Euler %%%% x = linspace(x0,xN,N); y = zeros(1,N); h = (xN-x0)/N; y(1) = 1; for k=2:N y(k) = y(k-1) + dydx(x(…

を考える。 clear; close all; Ns=[8,10,15,20,30]; x0=0; xN=20; figure(); hold on; for idx=1:length(Ns) N = Ns(idx); x=linspace(x0,xN,N); y=zeros(1,N); y(1)=0.1; h=(xN-x0)/N; f = @(y) y-y^2; for k=2:N y(k) = y(k-1) + h*(f(y(k-1))); end plot(…

クォータニオンの微分方程式 Matlabで解析解が出るっぽい >> syms w_x w_y w_z >> syms q0(t) q1(t) q2(t) q3(t) >> A=[0 -w_x -w_y -w_z;w_x 0 w_z -w_y; w_y -w_z 0 w_x; w_z w_y -w_x 0] A = [ 0, -w_x, -w_y, -w_z] [ w_x, 0, w_z, -w_y] [ w_y, -w_z, 0…