Why the results of three codes with the same expression and method are different ?
I used Newton-Raphson method to get velocity and temperature profiles.
I attached the differential equations as a file.
The desired result can only be obtained in the first code, and the specific matrix can be obtained in the remaining code, so the iteration value increases or the result itself cannot be obtained.
I will attach the ideal result as a graph.
code # 1
format longE
%% Predesignation
iter = 0; a = 0; b = 1; m = 1000; h = (b – a)/m; x = a: h : b; tol = 1e-7;
v0 = 1; vn = 40; t0 = 1;
ta = 0.543; k = 5.8;
st = [0.1, 0.5, 1]’;
max_iter = 1000;
for ii = 1 : 1 : length(st)
iter = 0;
iter(ii) = iter;
v_temp = v0:(40-1)/m:vn;
v = v_temp(2:end-1)’;
t_temp = 1:(0.9-1)/m:0.9;
t = t_temp(2:end)’;
v_t = [v ; t];
St = st(ii);
while iter <= max_iter
% Residual function
R1 = zeros(m – 1, 1);
for i = 1 : m – 1
if i == 1
R1(i) = 4*v(i)*v(i+1)*t(i)^2-8*v(i)^2*t(i)^2+4*1*v(i)*t(i)^2 …
-v(i+1)^2*t(i)^2+2*1*v(i+1)*t(i)^2 …
-1^2*t(i)^2-2*k*St*h*v(i)^(7/6)*(ta-t(i))*(v(i+1)-1);
elseif i == m – 1
R1(i) = 4*v(i)*40*t(i)^2-8*v(i)^2*t(i)^2+4*v(i-1)*v(i)*t(i)^2 …
-40^2*t(i)^2+2*v(i-1)*40*t(i)^2 …
-v(i-1)^2*t(i)^2-2*k*St*h*v(i)^(7/6)*(ta-t(i))*(40-v(i-1));
else
R1(i) = 4*v(i)*v(i+1)*t(i)^2-8*v(i)^2*t(i)^2+4*v(i-1)*v(i)*t(i)^2 …
-v(i+1)^2*t(i)^2+2*v(i-1)*v(i+1)*t(i)^2 …
-v(i-1)^2*t(i)^2-2*k*St*h*v(i)^(7/6)*(ta-t(i))*(v(i+1)-v(i-1));
end
end
R2 = zeros(m, 1);
for i = 1 : m
if i == 1
R2(i) = t(i+1)-1-2*h*St*v(i)^(1/6)*(ta-1);
elseif i == m
R2(i) = t(i-2)-4*t(i-1)+3*t(i)-2*h*St*40^(1/6)*(ta-t(i));
else
R2(i) = t(i+1)-t(i-1)-2*h*St*v(i)^(1/6)*(ta-t(i));
end
end
res = [R1 ; R2];
% Jacobian
J11 = zeros(m – 1);
for i = 1 : m – 1
if i == 1
J11(i, i) = 4*v(i+1)*t(i)^2-16*v(i)*t(i)^2+4*1*t(i)^2-2*k*St*h*(ta-t(i))*(v(i+1)-1)*(7/6*v(i)^(1/6));
J11(i, i + 1) = 4*v(i)*t(i)^2-2*v(i+1)*t(i)^2+2*1*t(i)^2-2*k*St*h*v(i)^(7/6)*(ta-t(i));
elseif i == m – 1
J11(i, i – 1) = 4*v(i)*t(i)^2+2*40*t(i)^2-2*v(i-1)*t(i)^2+2*k*St*h*v(i)^(7/6)*(ta-t(i));
J11(i, i) = 4*40*t(i)^2-16*v(i)*t(i)^2+4*v(i-1)*t(i)^2-2*k*St*h*(ta-t(i))*(40-v(i-1))*(7/6*v(i)^(1/6));
else
J11(i, i – 1) = 4*v(i)*t(i)^2+2*v(i+1)*t(i)^2-2*v(i-1)*t(i)^2+2*k*St*h*v(i)^(7/6)*(ta-t(i));
J11(i, i) = 4*v(i+1)*t(i)^2-16*v(i)*t(i)^2+4*v(i-1)*t(i)^2-2*k*St*h*(ta-t(i))*(v(i+1)-v(i-1))*(7/6*v(i)^(1/6));
J11(i, i + 1) = 4*v(i)*t(i)^2-2*v(i+1)*t(i)^2+2*v(i-1)*t(i)^2-2*k*St*h*v(i)^(7/6)*(ta-t(i));
end
end
J12 = zeros(m – 1, m);
for i = 1 : m – 1
if i == 1
J12(i, i) = 8*v(i)*v(i+1)*t(i)-16*v(i)^2*t(i)+8*1*v(i)*t(i)-2*v(i+1)^2*t(i) …
+4*1*v(i+1)*t(i)-2*1^2*t(i)+2*k*St*h*v(i)^(7/6)*(v(i+1)-1);
elseif i == m – 1
J12(i, i) = 8*v(i)*40*t(i)-16*v(i)^2*t(i)+8*v(i-1)*v(i)*t(i)-2*40^2*t(i) …
+4*v(i-1)*40*t(i)-2*v(i-1)^2*t(i)+2*k*St*h*v(i)^(7/6)*(40-v(i-1));
else
J12(i, i) = 8*v(i)*v(i+1)*t(i)-16*v(i)^2*t(i)+8*v(i-1)*v(i)*t(i)-2*v(i+1)^2*t(i) …
+4*v(i-1)*v(i+1)*t(i)-2*v(i-1)^2*t(i)+2*k*St*h*v(i)^(7/6)*(v(i+1)-v(i-1));
end
end
J21 = zeros(m, m – 1);
for i = 1 : m – 1
J21(i, i) = -1/3*h*St*(ta-t(i))*v(i)^(-5/6);
end
J22 = zeros(m, m);
for i = 1 : m
if i == 1
J22(i, i) = 2*h*St*v(i)^(1/6);
J22(i, i + 1) = 1;
elseif i == m
J22(i, i – 2) = 1;
J22(i, i – 1) = -4;
J22(i, i) = 3+2*h*St*40^(1/6);
else
J22(i, i – 1) = -1;
J22(i, i) = 2*h*St*v(i)^(1/6);
J22(i, i + 1) = 1;
end
end
jac = [J11, J12 ; J21, J22];
%% Updating
delv_delt = -inv(jac)*res; % res를 처음부터 column vector 형태로 만들자 다른점 하나 추가
error = sqrt(sum(delv_delt.^2));
%% Escaping condition
if error <= tol
v = v_t(1 : m – 1);
t = v_t(m : end);
v_compile(:, ii) = [v0 ; v ; vn];
t_compile(:, ii) = [t0 ; t];
break
else
iter = iter + 1;
v_t = v_t + delv_delt;
v = v_t(1 : m – 1);
t = v_t(m : end);
end
end
plot(x, v_compile(:, ii), ‘r’)
hold on
plot(x, t_compile(:, ii), ‘b’)
hold on
end
code #2
format longE
%% Predesignation
iter = 0; a = 0; b = 1; m = 1000; h = (b – a)/m; x = a: h : b; tol = 1e-7;
v0 = 1; vn = 40; t0 = 1;
ta = 0.543; k = 5.8;
st = [0.1, 0.5, 1]’;
max_iter = 1000;
v_temp = v0:(40-1)/m:vn;
v = v_temp(2:end-1)’;
t_temp = 1:(0.9-1)/m:0.9;
t = t_temp(2:end)’;
for ii = 1 : 1 : length(st)
iter = 0;
iter(ii) = iter;
v_t = [v ; t];
St = st(ii);
while iter <= max_iter
% Residual function
R1 = zeros(m – 1, 1);
for i = 1 : m – 1
if i == 1
R1(i) = 4*v(i)*v(i+1)*t(i)^2-8*v(i)^2*t(i)^2+4*1*v(i)*t(i)^2 …
-v(i+1)^2*t(i)^2+2*1*v(i+1)*t(i)^2 …
-1^2*t(i)^2-2*k*St*h*v(i)^(7/6)*(ta-t(i))*(v(i+1)-1);
elseif i == m – 1
R1(i) = 4*v(i)*40*t(i)^2-8*v(i)^2*t(i)^2+4*v(i-1)*v(i)*t(i)^2 …
-40^2*t(i)^2+2*v(i-1)*40*t(i)^2 …
-v(i-1)^2*t(i)^2-2*k*St*h*v(i)^(7/6)*(ta-t(i))*(40-v(i-1));
else
R1(i) = 4*v(i)*v(i+1)*t(i)^2-8*v(i)^2*t(i)^2+4*v(i-1)*v(i)*t(i)^2 …
-v(i+1)^2*t(i)^2+2*v(i-1)*v(i+1)*t(i)^2 …
-v(i-1)^2*t(i)^2-2*k*St*h*v(i)^(7/6)*(ta-t(i))*(v(i+1)-v(i-1));
end
end
R2 = zeros(m, 1);
for i = 1 : m
if i == 1
R2(i) = t(i+1)-1-2*h*St*v(i)^(1/6)*(ta-1);
elseif i == m
R2(i) = t(i-2)-4*t(i-1)+3*t(i)-2*h*St*40^(1/6)*(ta-t(i));
else
R2(i) = t(i+1)-t(i-1)-2*h*St*v(i)^(1/6)*(ta-t(i));
end
end
res = [R1 ; R2];
% Jacobian
J11 = zeros(m – 1);
for i = 1 : m – 1
if i == 1
J11(i, i) = 4*v(i+1)*t(i)^2-16*v(i)*t(i)^2+4*1*t(i)^2-2*k*St*h*(ta-t(i))*(v(i+1)-1)*(7/6*v(i)^(1/6));
J11(i, i + 1) = 4*v(i)*t(i)^2-2*v(i+1)*t(i)^2+2*1*t(i)^2-2*k*St*h*v(i)^(7/6)*(ta-t(i));
elseif i == m – 1
J11(i, i – 1) = 4*v(i)*t(i)^2+2*40*t(i)^2-2*v(i-1)*t(i)^2+2*k*St*h*v(i)^(7/6)*(ta-t(i));
J11(i, i) = 4*40*t(i)^2-16*v(i)*t(i)^2+4*v(i-1)*t(i)^2-2*k*St*h*(ta-t(i))*(40-v(i-1))*(7/6*v(i)^(1/6));
else
J11(i, i – 1) = 4*v(i)*t(i)^2+2*v(i+1)*t(i)^2-2*v(i-1)*t(i)^2+2*k*St*h*v(i)^(7/6)*(ta-t(i));
J11(i, i) = 4*v(i+1)*t(i)^2-16*v(i)*t(i)^2+4*v(i-1)*t(i)^2-2*k*St*h*(ta-t(i))*(v(i+1)-v(i-1))*(7/6*v(i)^(1/6));
J11(i, i + 1) = 4*v(i)*t(i)^2-2*v(i+1)*t(i)^2+2*v(i-1)*t(i)^2-2*k*St*h*v(i)^(7/6)*(ta-t(i));
end
end
J12 = zeros(m – 1, m);
for i = 1 : m – 1
if i == 1
J12(i, i) = 8*v(i)*v(i+1)*t(i)-16*v(i)^2*t(i)+8*1*v(i)*t(i)-2*v(i+1)^2*t(i) …
+4*1*v(i+1)*t(i)-2*1^2*t(i)+2*k*St*h*v(i)^(7/6)*(v(i+1)-1);
elseif i == m – 1
J12(i, i) = 8*v(i)*40*t(i)-16*v(i)^2*t(i)+8*v(i-1)*v(i)*t(i)-2*40^2*t(i) …
+4*v(i-1)*40*t(i)-2*v(i-1)^2*t(i)+2*k*St*h*v(i)^(7/6)*(40-v(i-1));
else
J12(i, i) = 8*v(i)*v(i+1)*t(i)-16*v(i)^2*t(i)+8*v(i-1)*v(i)*t(i)-2*v(i+1)^2*t(i) …
+4*v(i-1)*v(i+1)*t(i)-2*v(i-1)^2*t(i)+2*k*St*h*v(i)^(7/6)*(v(i+1)-v(i-1));
end
end
J21 = zeros(m, m – 1);
for i = 1 : m – 1
J21(i, i) = -1/3*h*St*(ta-t(i))*v(i)^(-5/6);
end
J22 = zeros(m, m);
for i = 1 : m
if i == 1
J22(i, i) = 2*h*St*v(i)^(1/6);
J22(i, i + 1) = 1;
elseif i == m
J22(i, i – 2) = 1;
J22(i, i – 1) = -4;
J22(i, i) = 3+2*h*St*40^(1/6);
else
J22(i, i – 1) = -1;
J22(i, i) = 2*h*St*v(i)^(1/6);
J22(i, i + 1) = 1;
end
end
jac = [J11, J12 ; J21, J22];
%% Updating
delv_delt = -inv(jac)*res; % res를 처음부터 column vector 형태로 만들자 다른점 하나 추가
error = sqrt(sum(delv_delt.^2));
%% Escaping condition
if error <= tol
v = v_t(1 : m – 1);
t = v_t(m : end);
v_compile(:, ii) = [v0 ; v ; vn];
t_compile(:, ii) = [t0 ; t];
break
else
iter = iter + 1;
v_t = v_t + delv_delt;
v = v_t(1 : m – 1);
t = v_t(m : end);
end
end
plot(x, v_compile(:, ii), ‘r’)
hold on
plot(x, t_compile(:, ii), ‘b’)
hold on
end
code #3
format longE
%% Predesignation
a = 0; b = 1; m = 1000; h = (b – a)/m; x = a : h : b; max_iter = 300; tol = 1e-2;
vi = 1; ve = 40; ti = 1;
k = 5.8; ta = 0.543; St_v = [0.1 ; 0.5 ; 1];
%% for loop about St and 3rd dim. of v, t
for p = 1 : 1 : length(St_v)
j = 1;
iter = 0;
iter(p) = iter;
St = St_v(p);
if p == 1
v(:, j, p) = (vi : (ve – vi)/m : ve)’;
t(:, j, p) = (ti : (-1 – ti)/m : -1)’;
end
%% while loop
while iter <= max_iter
%% Residual function
R1 = zeros(m – 1, 1);
R1(1, 1) = -2*h*k*v(2, j, p)^(7/6)*St*(ta – t(2, j, p))*(v(3, j, p) – v(1, j, p)) – t(2, j, p)^2*(v(3, j, p) – v(1, j, p))^2 + 4*t(2, j, p)^2*v(2, j, p)*(v(3, j, p) – 2*v(2, j, p) + v(1, j, p));
R1(end, 1) = -2*h*k*v(end – 1, j, p)^(7/6)*St*(ta – t(end – 1, j, p))*(v(end, j, p) – v(end – 2, j, p)) – t(end – 1, j, p)^2*(v(end, j, p) – v(end – 2, j, p))^2 + 4*t(end – 1, j, p)^2*v(end – 1, j, p)*(v(end, j, p) – 2*v(end – 1, j, p) + v(end – 2, j, p));
for i = 3 : 1 : m – 1
R1(i – 1, 1) = -2*h*k*v(i, j, p)^(7/6)*St*(ta – t(i, j, p))*(v(i + 1, j, p) – v(i – 1, j, p)) – t(i)^2*(v(i + 1, j, p) – v(i – 1, j, p))^2 + 4*t(i, j, p)^2*v(i, j, p)*(v(i + 1, j, p) – 2*v(i, j, p) + v(i – 1, j, p));
end
R2 = zeros(m, 1);
R2(1, 1) = t(3, j, p) – t(1, j, p) – 2*h*St*v(2, j, p)^(1/6)*(ta – t(2, j, p));
R2(end, 1) = 3*t(end, j, p) – 4*t(end – 1, j, p) + t(end – 2, j, p) – 2*h*St*v(end, j, p)^(1/6)*(ta – t(end, j, p));
for i = 3 : 1 : m
R2(i – 1, 1) = t(i + 1, j, p) – t(i – 1, j, p) – 2*h*St*v(i, j, p)^(1/6)*(ta – t(i, j, p));
end
R = [R1 ; R2];
%% Jacobian
J11 = zeros(m – 1, m – 1);
J11(1, 1) = -2*h*k*(7/6)*v(2, j, p)^(1/6)*St*(ta – t(2, j, p))*(v(3, j, p) – v(1, j, p)) + 4*t(2, j, p)^2*(v(3, j, p) – 2*v(2, j, p) + v(1, j, p)) + 4*t(2, j, p)^2*v(2, j, p)*-2;
J11(1, 2) = -2*h*k*v(2, j, p)^(7/6)*St*(ta – t(2, j, p)) – t(2, j, p)^2*(2*v(3, j, p) – 2*v(1, j, p)) + 4*t(2, j, p)^2*v(2, j, p);
J11(end, end – 1) = 2*h*k*v(end – 1, j , p)^(7/6)*St*(ta – t(end – 1, j, p)) – t(end – 1, j , p)^2*(2*v(end – 2, j, p) – 2*v(end, j, p)) + 4*t(end – 1, j, p)^2*v(end – 1, j, p);
J11(end, end) = -2*h*k*(7/6)*v(end – 1, j, p)^(1/6)*St*(ta – t(end – 1, j , p))*(v(end, j, p) – v(end – 2, j, p)) + 4*t(end – 1, j, p)^2*(v(end, j, p) – 2*v(end – 1, j, p) + v(end – 2, j, p)) + 4*t(end – 1, j, p)^2*v(end – 1, j, p)*-2;
for i = 3 : 1 : m – 1
J11(i – 1, i – 2) = 2*h*k*v(i, j, p)^(7/6)*St*(ta – t(i, j, p)) – t(i)^2*(2*v(i – 1, j, p) – 2*v(i + 1, j, p)) + 4*t(i, j, p)^2*v(i, j, p);
J11(i – 1, i – 1) = -2*h*k*(7/6)*v(i, j, p)^(1/6)*St*(ta – t(i, j, p))*(v(i + 1, j, p) – v(i – 1, j, p)) + 4*t(i, j, p)^2*(v(i + 1, j, p) – 2*v(i, j, p) + v(i – 1, j, p)) + 4*t(i, j, p)^2*v(i, j, p)*-2;
J11(i – 1, i) = -2*h*k*v(i, j, p)^(7/6)*St*(ta – t(i, j, p)) – t(i, j, p)^2*(2*v(i + 1, j, p) – 2*v(i – 1, j, p)) + 4*t(i, j, p)^2*v(i, j, p);
end
J12 = zeros(m – 1, m);
J12(1, 1) = 2*h*k*v(2, j , p)^(7/6)*St*(v(3, j, p) – v(1, j, p)) – 2*t(2, j, p)*(v(3, j, p) – v(1, j, p))^2 + 8*t(2, j, p)*v(2, j, p)*(v(3, j, p) – 2*v(2, j, p) + v(1, j, p));
J12(end, end – 1) = 2*h*k*v(end – 1, j, p)^(7/6)*St*(v(end, j, p) – v(end – 2, j, p)) – 2*t(end – 1, j, p)*(v(end, j, p) – v(end – 2, j, p))^2 + 8*t(end – 1, j, p)*v(end – 1, j, p)*(v(end, j, p) – 2*v(end – 1, j, p) + v(end – 2, j, p));
for i = 3 : 1 : m – 1
J12(i – 1, i – 1) = 2*h*k*v(i, j, p)^(7/6)*St*(v(i + 1, j, p) – v(i – 1, j, p)) – 2*t(i, j, p)*(v(i + 1, j, p) – v(i – 1, j, p))^2 + 8*t(i, j, p)*v(i, j, p)*(v(i + 1, j, p) – 2*v(i, j, p) + v(i – 1, j, p));
end
J21 = zeros(m, m – 1);
for i = 2 : 1 : m
J21(i – 1, i – 1) = -2*h*St*(1/6)*v(i, j, p)^(-5/6)*(ta – t(i, j, p));
end
J22 = zeros(m, m);
J22(1, 1) = 2*h*St*v(2, j, p)^(1/6);
J22(1, 2) = 1;
J22(end, end – 2) = 1;
J22(end, end – 1) = -4;
J22(end, end) = 3 + 2*h*St*v(end, j, p)^(1/6);
for i = 3 : 1 : m
J22(i – 1, i – 2) = -1;
J22(i – 1, i – 1) = 2*h*St*v(i, j, p)^(1/6);
J22(i – 1, i) = 1;
end
J = [J11 J12 ; J21 J22];
%% Updating & Condition number
delta = -inv(J)*R;
J_i = inv(J);
for i = 1 : 1 : length(J)
aa(i) = sum(abs(J(i, :)));
bb(i) = sum(abs((J_i(i, :))));
end
aam = max(aa);
bbm = max(bb);
c = aam*bbm;
C(j, p) = 1/(aam*bbm);
%% Escaping condition & Condition number of jacobian0
if sqrt(sum(delta.^2)) <= tol
vf(:, p) = v(:, j, p);
tf(:, p) = t(:, j, p);
if p <= length(St_v) – 1
j = 1;
v(:, j, p + 1) = vf(:, p);
t(:, j, p + 1) = tf(:, p);
else
v(:, j, p + 1) = 0;
t(:, j, p + 1) = 0;
end
break
else
v(:, j + 1, p) = v(:, j, p);
v(2 : 1 : end – 1, j + 1, p) = v(2 : 1 : end – 1, j, p) + delta(1 : 1 : m – 1);
t(:, j + 1, p) = t(:, j, p);
t(2 : 1 : end, j + 1, p) = t(2 : 1 : end, j, p) + delta(m : 1 : end);
j = j + 1;
iter = iter + 1;
end
end
%% Plotting
plot(x, vf(:, p), ‘r’)
hold on
plot(x, tf(:, p), ‘b’)
hold on
endI used Newton-Raphson method to get velocity and temperature profiles.
I attached the differential equations as a file.
The desired result can only be obtained in the first code, and the specific matrix can be obtained in the remaining code, so the iteration value increases or the result itself cannot be obtained.
I will attach the ideal result as a graph.
code # 1
format longE
%% Predesignation
iter = 0; a = 0; b = 1; m = 1000; h = (b – a)/m; x = a: h : b; tol = 1e-7;
v0 = 1; vn = 40; t0 = 1;
ta = 0.543; k = 5.8;
st = [0.1, 0.5, 1]’;
max_iter = 1000;
for ii = 1 : 1 : length(st)
iter = 0;
iter(ii) = iter;
v_temp = v0:(40-1)/m:vn;
v = v_temp(2:end-1)’;
t_temp = 1:(0.9-1)/m:0.9;
t = t_temp(2:end)’;
v_t = [v ; t];
St = st(ii);
while iter <= max_iter
% Residual function
R1 = zeros(m – 1, 1);
for i = 1 : m – 1
if i == 1
R1(i) = 4*v(i)*v(i+1)*t(i)^2-8*v(i)^2*t(i)^2+4*1*v(i)*t(i)^2 …
-v(i+1)^2*t(i)^2+2*1*v(i+1)*t(i)^2 …
-1^2*t(i)^2-2*k*St*h*v(i)^(7/6)*(ta-t(i))*(v(i+1)-1);
elseif i == m – 1
R1(i) = 4*v(i)*40*t(i)^2-8*v(i)^2*t(i)^2+4*v(i-1)*v(i)*t(i)^2 …
-40^2*t(i)^2+2*v(i-1)*40*t(i)^2 …
-v(i-1)^2*t(i)^2-2*k*St*h*v(i)^(7/6)*(ta-t(i))*(40-v(i-1));
else
R1(i) = 4*v(i)*v(i+1)*t(i)^2-8*v(i)^2*t(i)^2+4*v(i-1)*v(i)*t(i)^2 …
-v(i+1)^2*t(i)^2+2*v(i-1)*v(i+1)*t(i)^2 …
-v(i-1)^2*t(i)^2-2*k*St*h*v(i)^(7/6)*(ta-t(i))*(v(i+1)-v(i-1));
end
end
R2 = zeros(m, 1);
for i = 1 : m
if i == 1
R2(i) = t(i+1)-1-2*h*St*v(i)^(1/6)*(ta-1);
elseif i == m
R2(i) = t(i-2)-4*t(i-1)+3*t(i)-2*h*St*40^(1/6)*(ta-t(i));
else
R2(i) = t(i+1)-t(i-1)-2*h*St*v(i)^(1/6)*(ta-t(i));
end
end
res = [R1 ; R2];
% Jacobian
J11 = zeros(m – 1);
for i = 1 : m – 1
if i == 1
J11(i, i) = 4*v(i+1)*t(i)^2-16*v(i)*t(i)^2+4*1*t(i)^2-2*k*St*h*(ta-t(i))*(v(i+1)-1)*(7/6*v(i)^(1/6));
J11(i, i + 1) = 4*v(i)*t(i)^2-2*v(i+1)*t(i)^2+2*1*t(i)^2-2*k*St*h*v(i)^(7/6)*(ta-t(i));
elseif i == m – 1
J11(i, i – 1) = 4*v(i)*t(i)^2+2*40*t(i)^2-2*v(i-1)*t(i)^2+2*k*St*h*v(i)^(7/6)*(ta-t(i));
J11(i, i) = 4*40*t(i)^2-16*v(i)*t(i)^2+4*v(i-1)*t(i)^2-2*k*St*h*(ta-t(i))*(40-v(i-1))*(7/6*v(i)^(1/6));
else
J11(i, i – 1) = 4*v(i)*t(i)^2+2*v(i+1)*t(i)^2-2*v(i-1)*t(i)^2+2*k*St*h*v(i)^(7/6)*(ta-t(i));
J11(i, i) = 4*v(i+1)*t(i)^2-16*v(i)*t(i)^2+4*v(i-1)*t(i)^2-2*k*St*h*(ta-t(i))*(v(i+1)-v(i-1))*(7/6*v(i)^(1/6));
J11(i, i + 1) = 4*v(i)*t(i)^2-2*v(i+1)*t(i)^2+2*v(i-1)*t(i)^2-2*k*St*h*v(i)^(7/6)*(ta-t(i));
end
end
J12 = zeros(m – 1, m);
for i = 1 : m – 1
if i == 1
J12(i, i) = 8*v(i)*v(i+1)*t(i)-16*v(i)^2*t(i)+8*1*v(i)*t(i)-2*v(i+1)^2*t(i) …
+4*1*v(i+1)*t(i)-2*1^2*t(i)+2*k*St*h*v(i)^(7/6)*(v(i+1)-1);
elseif i == m – 1
J12(i, i) = 8*v(i)*40*t(i)-16*v(i)^2*t(i)+8*v(i-1)*v(i)*t(i)-2*40^2*t(i) …
+4*v(i-1)*40*t(i)-2*v(i-1)^2*t(i)+2*k*St*h*v(i)^(7/6)*(40-v(i-1));
else
J12(i, i) = 8*v(i)*v(i+1)*t(i)-16*v(i)^2*t(i)+8*v(i-1)*v(i)*t(i)-2*v(i+1)^2*t(i) …
+4*v(i-1)*v(i+1)*t(i)-2*v(i-1)^2*t(i)+2*k*St*h*v(i)^(7/6)*(v(i+1)-v(i-1));
end
end
J21 = zeros(m, m – 1);
for i = 1 : m – 1
J21(i, i) = -1/3*h*St*(ta-t(i))*v(i)^(-5/6);
end
J22 = zeros(m, m);
for i = 1 : m
if i == 1
J22(i, i) = 2*h*St*v(i)^(1/6);
J22(i, i + 1) = 1;
elseif i == m
J22(i, i – 2) = 1;
J22(i, i – 1) = -4;
J22(i, i) = 3+2*h*St*40^(1/6);
else
J22(i, i – 1) = -1;
J22(i, i) = 2*h*St*v(i)^(1/6);
J22(i, i + 1) = 1;
end
end
jac = [J11, J12 ; J21, J22];
%% Updating
delv_delt = -inv(jac)*res; % res를 처음부터 column vector 형태로 만들자 다른점 하나 추가
error = sqrt(sum(delv_delt.^2));
%% Escaping condition
if error <= tol
v = v_t(1 : m – 1);
t = v_t(m : end);
v_compile(:, ii) = [v0 ; v ; vn];
t_compile(:, ii) = [t0 ; t];
break
else
iter = iter + 1;
v_t = v_t + delv_delt;
v = v_t(1 : m – 1);
t = v_t(m : end);
end
end
plot(x, v_compile(:, ii), ‘r’)
hold on
plot(x, t_compile(:, ii), ‘b’)
hold on
end
code #2
format longE
%% Predesignation
iter = 0; a = 0; b = 1; m = 1000; h = (b – a)/m; x = a: h : b; tol = 1e-7;
v0 = 1; vn = 40; t0 = 1;
ta = 0.543; k = 5.8;
st = [0.1, 0.5, 1]’;
max_iter = 1000;
v_temp = v0:(40-1)/m:vn;
v = v_temp(2:end-1)’;
t_temp = 1:(0.9-1)/m:0.9;
t = t_temp(2:end)’;
for ii = 1 : 1 : length(st)
iter = 0;
iter(ii) = iter;
v_t = [v ; t];
St = st(ii);
while iter <= max_iter
% Residual function
R1 = zeros(m – 1, 1);
for i = 1 : m – 1
if i == 1
R1(i) = 4*v(i)*v(i+1)*t(i)^2-8*v(i)^2*t(i)^2+4*1*v(i)*t(i)^2 …
-v(i+1)^2*t(i)^2+2*1*v(i+1)*t(i)^2 …
-1^2*t(i)^2-2*k*St*h*v(i)^(7/6)*(ta-t(i))*(v(i+1)-1);
elseif i == m – 1
R1(i) = 4*v(i)*40*t(i)^2-8*v(i)^2*t(i)^2+4*v(i-1)*v(i)*t(i)^2 …
-40^2*t(i)^2+2*v(i-1)*40*t(i)^2 …
-v(i-1)^2*t(i)^2-2*k*St*h*v(i)^(7/6)*(ta-t(i))*(40-v(i-1));
else
R1(i) = 4*v(i)*v(i+1)*t(i)^2-8*v(i)^2*t(i)^2+4*v(i-1)*v(i)*t(i)^2 …
-v(i+1)^2*t(i)^2+2*v(i-1)*v(i+1)*t(i)^2 …
-v(i-1)^2*t(i)^2-2*k*St*h*v(i)^(7/6)*(ta-t(i))*(v(i+1)-v(i-1));
end
end
R2 = zeros(m, 1);
for i = 1 : m
if i == 1
R2(i) = t(i+1)-1-2*h*St*v(i)^(1/6)*(ta-1);
elseif i == m
R2(i) = t(i-2)-4*t(i-1)+3*t(i)-2*h*St*40^(1/6)*(ta-t(i));
else
R2(i) = t(i+1)-t(i-1)-2*h*St*v(i)^(1/6)*(ta-t(i));
end
end
res = [R1 ; R2];
% Jacobian
J11 = zeros(m – 1);
for i = 1 : m – 1
if i == 1
J11(i, i) = 4*v(i+1)*t(i)^2-16*v(i)*t(i)^2+4*1*t(i)^2-2*k*St*h*(ta-t(i))*(v(i+1)-1)*(7/6*v(i)^(1/6));
J11(i, i + 1) = 4*v(i)*t(i)^2-2*v(i+1)*t(i)^2+2*1*t(i)^2-2*k*St*h*v(i)^(7/6)*(ta-t(i));
elseif i == m – 1
J11(i, i – 1) = 4*v(i)*t(i)^2+2*40*t(i)^2-2*v(i-1)*t(i)^2+2*k*St*h*v(i)^(7/6)*(ta-t(i));
J11(i, i) = 4*40*t(i)^2-16*v(i)*t(i)^2+4*v(i-1)*t(i)^2-2*k*St*h*(ta-t(i))*(40-v(i-1))*(7/6*v(i)^(1/6));
else
J11(i, i – 1) = 4*v(i)*t(i)^2+2*v(i+1)*t(i)^2-2*v(i-1)*t(i)^2+2*k*St*h*v(i)^(7/6)*(ta-t(i));
J11(i, i) = 4*v(i+1)*t(i)^2-16*v(i)*t(i)^2+4*v(i-1)*t(i)^2-2*k*St*h*(ta-t(i))*(v(i+1)-v(i-1))*(7/6*v(i)^(1/6));
J11(i, i + 1) = 4*v(i)*t(i)^2-2*v(i+1)*t(i)^2+2*v(i-1)*t(i)^2-2*k*St*h*v(i)^(7/6)*(ta-t(i));
end
end
J12 = zeros(m – 1, m);
for i = 1 : m – 1
if i == 1
J12(i, i) = 8*v(i)*v(i+1)*t(i)-16*v(i)^2*t(i)+8*1*v(i)*t(i)-2*v(i+1)^2*t(i) …
+4*1*v(i+1)*t(i)-2*1^2*t(i)+2*k*St*h*v(i)^(7/6)*(v(i+1)-1);
elseif i == m – 1
J12(i, i) = 8*v(i)*40*t(i)-16*v(i)^2*t(i)+8*v(i-1)*v(i)*t(i)-2*40^2*t(i) …
+4*v(i-1)*40*t(i)-2*v(i-1)^2*t(i)+2*k*St*h*v(i)^(7/6)*(40-v(i-1));
else
J12(i, i) = 8*v(i)*v(i+1)*t(i)-16*v(i)^2*t(i)+8*v(i-1)*v(i)*t(i)-2*v(i+1)^2*t(i) …
+4*v(i-1)*v(i+1)*t(i)-2*v(i-1)^2*t(i)+2*k*St*h*v(i)^(7/6)*(v(i+1)-v(i-1));
end
end
J21 = zeros(m, m – 1);
for i = 1 : m – 1
J21(i, i) = -1/3*h*St*(ta-t(i))*v(i)^(-5/6);
end
J22 = zeros(m, m);
for i = 1 : m
if i == 1
J22(i, i) = 2*h*St*v(i)^(1/6);
J22(i, i + 1) = 1;
elseif i == m
J22(i, i – 2) = 1;
J22(i, i – 1) = -4;
J22(i, i) = 3+2*h*St*40^(1/6);
else
J22(i, i – 1) = -1;
J22(i, i) = 2*h*St*v(i)^(1/6);
J22(i, i + 1) = 1;
end
end
jac = [J11, J12 ; J21, J22];
%% Updating
delv_delt = -inv(jac)*res; % res를 처음부터 column vector 형태로 만들자 다른점 하나 추가
error = sqrt(sum(delv_delt.^2));
%% Escaping condition
if error <= tol
v = v_t(1 : m – 1);
t = v_t(m : end);
v_compile(:, ii) = [v0 ; v ; vn];
t_compile(:, ii) = [t0 ; t];
break
else
iter = iter + 1;
v_t = v_t + delv_delt;
v = v_t(1 : m – 1);
t = v_t(m : end);
end
end
plot(x, v_compile(:, ii), ‘r’)
hold on
plot(x, t_compile(:, ii), ‘b’)
hold on
end
code #3
format longE
%% Predesignation
a = 0; b = 1; m = 1000; h = (b – a)/m; x = a : h : b; max_iter = 300; tol = 1e-2;
vi = 1; ve = 40; ti = 1;
k = 5.8; ta = 0.543; St_v = [0.1 ; 0.5 ; 1];
%% for loop about St and 3rd dim. of v, t
for p = 1 : 1 : length(St_v)
j = 1;
iter = 0;
iter(p) = iter;
St = St_v(p);
if p == 1
v(:, j, p) = (vi : (ve – vi)/m : ve)’;
t(:, j, p) = (ti : (-1 – ti)/m : -1)’;
end
%% while loop
while iter <= max_iter
%% Residual function
R1 = zeros(m – 1, 1);
R1(1, 1) = -2*h*k*v(2, j, p)^(7/6)*St*(ta – t(2, j, p))*(v(3, j, p) – v(1, j, p)) – t(2, j, p)^2*(v(3, j, p) – v(1, j, p))^2 + 4*t(2, j, p)^2*v(2, j, p)*(v(3, j, p) – 2*v(2, j, p) + v(1, j, p));
R1(end, 1) = -2*h*k*v(end – 1, j, p)^(7/6)*St*(ta – t(end – 1, j, p))*(v(end, j, p) – v(end – 2, j, p)) – t(end – 1, j, p)^2*(v(end, j, p) – v(end – 2, j, p))^2 + 4*t(end – 1, j, p)^2*v(end – 1, j, p)*(v(end, j, p) – 2*v(end – 1, j, p) + v(end – 2, j, p));
for i = 3 : 1 : m – 1
R1(i – 1, 1) = -2*h*k*v(i, j, p)^(7/6)*St*(ta – t(i, j, p))*(v(i + 1, j, p) – v(i – 1, j, p)) – t(i)^2*(v(i + 1, j, p) – v(i – 1, j, p))^2 + 4*t(i, j, p)^2*v(i, j, p)*(v(i + 1, j, p) – 2*v(i, j, p) + v(i – 1, j, p));
end
R2 = zeros(m, 1);
R2(1, 1) = t(3, j, p) – t(1, j, p) – 2*h*St*v(2, j, p)^(1/6)*(ta – t(2, j, p));
R2(end, 1) = 3*t(end, j, p) – 4*t(end – 1, j, p) + t(end – 2, j, p) – 2*h*St*v(end, j, p)^(1/6)*(ta – t(end, j, p));
for i = 3 : 1 : m
R2(i – 1, 1) = t(i + 1, j, p) – t(i – 1, j, p) – 2*h*St*v(i, j, p)^(1/6)*(ta – t(i, j, p));
end
R = [R1 ; R2];
%% Jacobian
J11 = zeros(m – 1, m – 1);
J11(1, 1) = -2*h*k*(7/6)*v(2, j, p)^(1/6)*St*(ta – t(2, j, p))*(v(3, j, p) – v(1, j, p)) + 4*t(2, j, p)^2*(v(3, j, p) – 2*v(2, j, p) + v(1, j, p)) + 4*t(2, j, p)^2*v(2, j, p)*-2;
J11(1, 2) = -2*h*k*v(2, j, p)^(7/6)*St*(ta – t(2, j, p)) – t(2, j, p)^2*(2*v(3, j, p) – 2*v(1, j, p)) + 4*t(2, j, p)^2*v(2, j, p);
J11(end, end – 1) = 2*h*k*v(end – 1, j , p)^(7/6)*St*(ta – t(end – 1, j, p)) – t(end – 1, j , p)^2*(2*v(end – 2, j, p) – 2*v(end, j, p)) + 4*t(end – 1, j, p)^2*v(end – 1, j, p);
J11(end, end) = -2*h*k*(7/6)*v(end – 1, j, p)^(1/6)*St*(ta – t(end – 1, j , p))*(v(end, j, p) – v(end – 2, j, p)) + 4*t(end – 1, j, p)^2*(v(end, j, p) – 2*v(end – 1, j, p) + v(end – 2, j, p)) + 4*t(end – 1, j, p)^2*v(end – 1, j, p)*-2;
for i = 3 : 1 : m – 1
J11(i – 1, i – 2) = 2*h*k*v(i, j, p)^(7/6)*St*(ta – t(i, j, p)) – t(i)^2*(2*v(i – 1, j, p) – 2*v(i + 1, j, p)) + 4*t(i, j, p)^2*v(i, j, p);
J11(i – 1, i – 1) = -2*h*k*(7/6)*v(i, j, p)^(1/6)*St*(ta – t(i, j, p))*(v(i + 1, j, p) – v(i – 1, j, p)) + 4*t(i, j, p)^2*(v(i + 1, j, p) – 2*v(i, j, p) + v(i – 1, j, p)) + 4*t(i, j, p)^2*v(i, j, p)*-2;
J11(i – 1, i) = -2*h*k*v(i, j, p)^(7/6)*St*(ta – t(i, j, p)) – t(i, j, p)^2*(2*v(i + 1, j, p) – 2*v(i – 1, j, p)) + 4*t(i, j, p)^2*v(i, j, p);
end
J12 = zeros(m – 1, m);
J12(1, 1) = 2*h*k*v(2, j , p)^(7/6)*St*(v(3, j, p) – v(1, j, p)) – 2*t(2, j, p)*(v(3, j, p) – v(1, j, p))^2 + 8*t(2, j, p)*v(2, j, p)*(v(3, j, p) – 2*v(2, j, p) + v(1, j, p));
J12(end, end – 1) = 2*h*k*v(end – 1, j, p)^(7/6)*St*(v(end, j, p) – v(end – 2, j, p)) – 2*t(end – 1, j, p)*(v(end, j, p) – v(end – 2, j, p))^2 + 8*t(end – 1, j, p)*v(end – 1, j, p)*(v(end, j, p) – 2*v(end – 1, j, p) + v(end – 2, j, p));
for i = 3 : 1 : m – 1
J12(i – 1, i – 1) = 2*h*k*v(i, j, p)^(7/6)*St*(v(i + 1, j, p) – v(i – 1, j, p)) – 2*t(i, j, p)*(v(i + 1, j, p) – v(i – 1, j, p))^2 + 8*t(i, j, p)*v(i, j, p)*(v(i + 1, j, p) – 2*v(i, j, p) + v(i – 1, j, p));
end
J21 = zeros(m, m – 1);
for i = 2 : 1 : m
J21(i – 1, i – 1) = -2*h*St*(1/6)*v(i, j, p)^(-5/6)*(ta – t(i, j, p));
end
J22 = zeros(m, m);
J22(1, 1) = 2*h*St*v(2, j, p)^(1/6);
J22(1, 2) = 1;
J22(end, end – 2) = 1;
J22(end, end – 1) = -4;
J22(end, end) = 3 + 2*h*St*v(end, j, p)^(1/6);
for i = 3 : 1 : m
J22(i – 1, i – 2) = -1;
J22(i – 1, i – 1) = 2*h*St*v(i, j, p)^(1/6);
J22(i – 1, i) = 1;
end
J = [J11 J12 ; J21 J22];
%% Updating & Condition number
delta = -inv(J)*R;
J_i = inv(J);
for i = 1 : 1 : length(J)
aa(i) = sum(abs(J(i, :)));
bb(i) = sum(abs((J_i(i, :))));
end
aam = max(aa);
bbm = max(bb);
c = aam*bbm;
C(j, p) = 1/(aam*bbm);
%% Escaping condition & Condition number of jacobian0
if sqrt(sum(delta.^2)) <= tol
vf(:, p) = v(:, j, p);
tf(:, p) = t(:, j, p);
if p <= length(St_v) – 1
j = 1;
v(:, j, p + 1) = vf(:, p);
t(:, j, p + 1) = tf(:, p);
else
v(:, j, p + 1) = 0;
t(:, j, p + 1) = 0;
end
break
else
v(:, j + 1, p) = v(:, j, p);
v(2 : 1 : end – 1, j + 1, p) = v(2 : 1 : end – 1, j, p) + delta(1 : 1 : m – 1);
t(:, j + 1, p) = t(:, j, p);
t(2 : 1 : end, j + 1, p) = t(2 : 1 : end, j, p) + delta(m : 1 : end);
j = j + 1;
iter = iter + 1;
end
end
%% Plotting
plot(x, vf(:, p), ‘r’)
hold on
plot(x, tf(:, p), ‘b’)
hold on
end I used Newton-Raphson method to get velocity and temperature profiles.
I attached the differential equations as a file.
The desired result can only be obtained in the first code, and the specific matrix can be obtained in the remaining code, so the iteration value increases or the result itself cannot be obtained.
I will attach the ideal result as a graph.
code # 1
format longE
%% Predesignation
iter = 0; a = 0; b = 1; m = 1000; h = (b – a)/m; x = a: h : b; tol = 1e-7;
v0 = 1; vn = 40; t0 = 1;
ta = 0.543; k = 5.8;
st = [0.1, 0.5, 1]’;
max_iter = 1000;
for ii = 1 : 1 : length(st)
iter = 0;
iter(ii) = iter;
v_temp = v0:(40-1)/m:vn;
v = v_temp(2:end-1)’;
t_temp = 1:(0.9-1)/m:0.9;
t = t_temp(2:end)’;
v_t = [v ; t];
St = st(ii);
while iter <= max_iter
% Residual function
R1 = zeros(m – 1, 1);
for i = 1 : m – 1
if i == 1
R1(i) = 4*v(i)*v(i+1)*t(i)^2-8*v(i)^2*t(i)^2+4*1*v(i)*t(i)^2 …
-v(i+1)^2*t(i)^2+2*1*v(i+1)*t(i)^2 …
-1^2*t(i)^2-2*k*St*h*v(i)^(7/6)*(ta-t(i))*(v(i+1)-1);
elseif i == m – 1
R1(i) = 4*v(i)*40*t(i)^2-8*v(i)^2*t(i)^2+4*v(i-1)*v(i)*t(i)^2 …
-40^2*t(i)^2+2*v(i-1)*40*t(i)^2 …
-v(i-1)^2*t(i)^2-2*k*St*h*v(i)^(7/6)*(ta-t(i))*(40-v(i-1));
else
R1(i) = 4*v(i)*v(i+1)*t(i)^2-8*v(i)^2*t(i)^2+4*v(i-1)*v(i)*t(i)^2 …
-v(i+1)^2*t(i)^2+2*v(i-1)*v(i+1)*t(i)^2 …
-v(i-1)^2*t(i)^2-2*k*St*h*v(i)^(7/6)*(ta-t(i))*(v(i+1)-v(i-1));
end
end
R2 = zeros(m, 1);
for i = 1 : m
if i == 1
R2(i) = t(i+1)-1-2*h*St*v(i)^(1/6)*(ta-1);
elseif i == m
R2(i) = t(i-2)-4*t(i-1)+3*t(i)-2*h*St*40^(1/6)*(ta-t(i));
else
R2(i) = t(i+1)-t(i-1)-2*h*St*v(i)^(1/6)*(ta-t(i));
end
end
res = [R1 ; R2];
% Jacobian
J11 = zeros(m – 1);
for i = 1 : m – 1
if i == 1
J11(i, i) = 4*v(i+1)*t(i)^2-16*v(i)*t(i)^2+4*1*t(i)^2-2*k*St*h*(ta-t(i))*(v(i+1)-1)*(7/6*v(i)^(1/6));
J11(i, i + 1) = 4*v(i)*t(i)^2-2*v(i+1)*t(i)^2+2*1*t(i)^2-2*k*St*h*v(i)^(7/6)*(ta-t(i));
elseif i == m – 1
J11(i, i – 1) = 4*v(i)*t(i)^2+2*40*t(i)^2-2*v(i-1)*t(i)^2+2*k*St*h*v(i)^(7/6)*(ta-t(i));
J11(i, i) = 4*40*t(i)^2-16*v(i)*t(i)^2+4*v(i-1)*t(i)^2-2*k*St*h*(ta-t(i))*(40-v(i-1))*(7/6*v(i)^(1/6));
else
J11(i, i – 1) = 4*v(i)*t(i)^2+2*v(i+1)*t(i)^2-2*v(i-1)*t(i)^2+2*k*St*h*v(i)^(7/6)*(ta-t(i));
J11(i, i) = 4*v(i+1)*t(i)^2-16*v(i)*t(i)^2+4*v(i-1)*t(i)^2-2*k*St*h*(ta-t(i))*(v(i+1)-v(i-1))*(7/6*v(i)^(1/6));
J11(i, i + 1) = 4*v(i)*t(i)^2-2*v(i+1)*t(i)^2+2*v(i-1)*t(i)^2-2*k*St*h*v(i)^(7/6)*(ta-t(i));
end
end
J12 = zeros(m – 1, m);
for i = 1 : m – 1
if i == 1
J12(i, i) = 8*v(i)*v(i+1)*t(i)-16*v(i)^2*t(i)+8*1*v(i)*t(i)-2*v(i+1)^2*t(i) …
+4*1*v(i+1)*t(i)-2*1^2*t(i)+2*k*St*h*v(i)^(7/6)*(v(i+1)-1);
elseif i == m – 1
J12(i, i) = 8*v(i)*40*t(i)-16*v(i)^2*t(i)+8*v(i-1)*v(i)*t(i)-2*40^2*t(i) …
+4*v(i-1)*40*t(i)-2*v(i-1)^2*t(i)+2*k*St*h*v(i)^(7/6)*(40-v(i-1));
else
J12(i, i) = 8*v(i)*v(i+1)*t(i)-16*v(i)^2*t(i)+8*v(i-1)*v(i)*t(i)-2*v(i+1)^2*t(i) …
+4*v(i-1)*v(i+1)*t(i)-2*v(i-1)^2*t(i)+2*k*St*h*v(i)^(7/6)*(v(i+1)-v(i-1));
end
end
J21 = zeros(m, m – 1);
for i = 1 : m – 1
J21(i, i) = -1/3*h*St*(ta-t(i))*v(i)^(-5/6);
end
J22 = zeros(m, m);
for i = 1 : m
if i == 1
J22(i, i) = 2*h*St*v(i)^(1/6);
J22(i, i + 1) = 1;
elseif i == m
J22(i, i – 2) = 1;
J22(i, i – 1) = -4;
J22(i, i) = 3+2*h*St*40^(1/6);
else
J22(i, i – 1) = -1;
J22(i, i) = 2*h*St*v(i)^(1/6);
J22(i, i + 1) = 1;
end
end
jac = [J11, J12 ; J21, J22];
%% Updating
delv_delt = -inv(jac)*res; % res를 처음부터 column vector 형태로 만들자 다른점 하나 추가
error = sqrt(sum(delv_delt.^2));
%% Escaping condition
if error <= tol
v = v_t(1 : m – 1);
t = v_t(m : end);
v_compile(:, ii) = [v0 ; v ; vn];
t_compile(:, ii) = [t0 ; t];
break
else
iter = iter + 1;
v_t = v_t + delv_delt;
v = v_t(1 : m – 1);
t = v_t(m : end);
end
end
plot(x, v_compile(:, ii), ‘r’)
hold on
plot(x, t_compile(:, ii), ‘b’)
hold on
end
code #2
format longE
%% Predesignation
iter = 0; a = 0; b = 1; m = 1000; h = (b – a)/m; x = a: h : b; tol = 1e-7;
v0 = 1; vn = 40; t0 = 1;
ta = 0.543; k = 5.8;
st = [0.1, 0.5, 1]’;
max_iter = 1000;
v_temp = v0:(40-1)/m:vn;
v = v_temp(2:end-1)’;
t_temp = 1:(0.9-1)/m:0.9;
t = t_temp(2:end)’;
for ii = 1 : 1 : length(st)
iter = 0;
iter(ii) = iter;
v_t = [v ; t];
St = st(ii);
while iter <= max_iter
% Residual function
R1 = zeros(m – 1, 1);
for i = 1 : m – 1
if i == 1
R1(i) = 4*v(i)*v(i+1)*t(i)^2-8*v(i)^2*t(i)^2+4*1*v(i)*t(i)^2 …
-v(i+1)^2*t(i)^2+2*1*v(i+1)*t(i)^2 …
-1^2*t(i)^2-2*k*St*h*v(i)^(7/6)*(ta-t(i))*(v(i+1)-1);
elseif i == m – 1
R1(i) = 4*v(i)*40*t(i)^2-8*v(i)^2*t(i)^2+4*v(i-1)*v(i)*t(i)^2 …
-40^2*t(i)^2+2*v(i-1)*40*t(i)^2 …
-v(i-1)^2*t(i)^2-2*k*St*h*v(i)^(7/6)*(ta-t(i))*(40-v(i-1));
else
R1(i) = 4*v(i)*v(i+1)*t(i)^2-8*v(i)^2*t(i)^2+4*v(i-1)*v(i)*t(i)^2 …
-v(i+1)^2*t(i)^2+2*v(i-1)*v(i+1)*t(i)^2 …
-v(i-1)^2*t(i)^2-2*k*St*h*v(i)^(7/6)*(ta-t(i))*(v(i+1)-v(i-1));
end
end
R2 = zeros(m, 1);
for i = 1 : m
if i == 1
R2(i) = t(i+1)-1-2*h*St*v(i)^(1/6)*(ta-1);
elseif i == m
R2(i) = t(i-2)-4*t(i-1)+3*t(i)-2*h*St*40^(1/6)*(ta-t(i));
else
R2(i) = t(i+1)-t(i-1)-2*h*St*v(i)^(1/6)*(ta-t(i));
end
end
res = [R1 ; R2];
% Jacobian
J11 = zeros(m – 1);
for i = 1 : m – 1
if i == 1
J11(i, i) = 4*v(i+1)*t(i)^2-16*v(i)*t(i)^2+4*1*t(i)^2-2*k*St*h*(ta-t(i))*(v(i+1)-1)*(7/6*v(i)^(1/6));
J11(i, i + 1) = 4*v(i)*t(i)^2-2*v(i+1)*t(i)^2+2*1*t(i)^2-2*k*St*h*v(i)^(7/6)*(ta-t(i));
elseif i == m – 1
J11(i, i – 1) = 4*v(i)*t(i)^2+2*40*t(i)^2-2*v(i-1)*t(i)^2+2*k*St*h*v(i)^(7/6)*(ta-t(i));
J11(i, i) = 4*40*t(i)^2-16*v(i)*t(i)^2+4*v(i-1)*t(i)^2-2*k*St*h*(ta-t(i))*(40-v(i-1))*(7/6*v(i)^(1/6));
else
J11(i, i – 1) = 4*v(i)*t(i)^2+2*v(i+1)*t(i)^2-2*v(i-1)*t(i)^2+2*k*St*h*v(i)^(7/6)*(ta-t(i));
J11(i, i) = 4*v(i+1)*t(i)^2-16*v(i)*t(i)^2+4*v(i-1)*t(i)^2-2*k*St*h*(ta-t(i))*(v(i+1)-v(i-1))*(7/6*v(i)^(1/6));
J11(i, i + 1) = 4*v(i)*t(i)^2-2*v(i+1)*t(i)^2+2*v(i-1)*t(i)^2-2*k*St*h*v(i)^(7/6)*(ta-t(i));
end
end
J12 = zeros(m – 1, m);
for i = 1 : m – 1
if i == 1
J12(i, i) = 8*v(i)*v(i+1)*t(i)-16*v(i)^2*t(i)+8*1*v(i)*t(i)-2*v(i+1)^2*t(i) …
+4*1*v(i+1)*t(i)-2*1^2*t(i)+2*k*St*h*v(i)^(7/6)*(v(i+1)-1);
elseif i == m – 1
J12(i, i) = 8*v(i)*40*t(i)-16*v(i)^2*t(i)+8*v(i-1)*v(i)*t(i)-2*40^2*t(i) …
+4*v(i-1)*40*t(i)-2*v(i-1)^2*t(i)+2*k*St*h*v(i)^(7/6)*(40-v(i-1));
else
J12(i, i) = 8*v(i)*v(i+1)*t(i)-16*v(i)^2*t(i)+8*v(i-1)*v(i)*t(i)-2*v(i+1)^2*t(i) …
+4*v(i-1)*v(i+1)*t(i)-2*v(i-1)^2*t(i)+2*k*St*h*v(i)^(7/6)*(v(i+1)-v(i-1));
end
end
J21 = zeros(m, m – 1);
for i = 1 : m – 1
J21(i, i) = -1/3*h*St*(ta-t(i))*v(i)^(-5/6);
end
J22 = zeros(m, m);
for i = 1 : m
if i == 1
J22(i, i) = 2*h*St*v(i)^(1/6);
J22(i, i + 1) = 1;
elseif i == m
J22(i, i – 2) = 1;
J22(i, i – 1) = -4;
J22(i, i) = 3+2*h*St*40^(1/6);
else
J22(i, i – 1) = -1;
J22(i, i) = 2*h*St*v(i)^(1/6);
J22(i, i + 1) = 1;
end
end
jac = [J11, J12 ; J21, J22];
%% Updating
delv_delt = -inv(jac)*res; % res를 처음부터 column vector 형태로 만들자 다른점 하나 추가
error = sqrt(sum(delv_delt.^2));
%% Escaping condition
if error <= tol
v = v_t(1 : m – 1);
t = v_t(m : end);
v_compile(:, ii) = [v0 ; v ; vn];
t_compile(:, ii) = [t0 ; t];
break
else
iter = iter + 1;
v_t = v_t + delv_delt;
v = v_t(1 : m – 1);
t = v_t(m : end);
end
end
plot(x, v_compile(:, ii), ‘r’)
hold on
plot(x, t_compile(:, ii), ‘b’)
hold on
end
code #3
format longE
%% Predesignation
a = 0; b = 1; m = 1000; h = (b – a)/m; x = a : h : b; max_iter = 300; tol = 1e-2;
vi = 1; ve = 40; ti = 1;
k = 5.8; ta = 0.543; St_v = [0.1 ; 0.5 ; 1];
%% for loop about St and 3rd dim. of v, t
for p = 1 : 1 : length(St_v)
j = 1;
iter = 0;
iter(p) = iter;
St = St_v(p);
if p == 1
v(:, j, p) = (vi : (ve – vi)/m : ve)’;
t(:, j, p) = (ti : (-1 – ti)/m : -1)’;
end
%% while loop
while iter <= max_iter
%% Residual function
R1 = zeros(m – 1, 1);
R1(1, 1) = -2*h*k*v(2, j, p)^(7/6)*St*(ta – t(2, j, p))*(v(3, j, p) – v(1, j, p)) – t(2, j, p)^2*(v(3, j, p) – v(1, j, p))^2 + 4*t(2, j, p)^2*v(2, j, p)*(v(3, j, p) – 2*v(2, j, p) + v(1, j, p));
R1(end, 1) = -2*h*k*v(end – 1, j, p)^(7/6)*St*(ta – t(end – 1, j, p))*(v(end, j, p) – v(end – 2, j, p)) – t(end – 1, j, p)^2*(v(end, j, p) – v(end – 2, j, p))^2 + 4*t(end – 1, j, p)^2*v(end – 1, j, p)*(v(end, j, p) – 2*v(end – 1, j, p) + v(end – 2, j, p));
for i = 3 : 1 : m – 1
R1(i – 1, 1) = -2*h*k*v(i, j, p)^(7/6)*St*(ta – t(i, j, p))*(v(i + 1, j, p) – v(i – 1, j, p)) – t(i)^2*(v(i + 1, j, p) – v(i – 1, j, p))^2 + 4*t(i, j, p)^2*v(i, j, p)*(v(i + 1, j, p) – 2*v(i, j, p) + v(i – 1, j, p));
end
R2 = zeros(m, 1);
R2(1, 1) = t(3, j, p) – t(1, j, p) – 2*h*St*v(2, j, p)^(1/6)*(ta – t(2, j, p));
R2(end, 1) = 3*t(end, j, p) – 4*t(end – 1, j, p) + t(end – 2, j, p) – 2*h*St*v(end, j, p)^(1/6)*(ta – t(end, j, p));
for i = 3 : 1 : m
R2(i – 1, 1) = t(i + 1, j, p) – t(i – 1, j, p) – 2*h*St*v(i, j, p)^(1/6)*(ta – t(i, j, p));
end
R = [R1 ; R2];
%% Jacobian
J11 = zeros(m – 1, m – 1);
J11(1, 1) = -2*h*k*(7/6)*v(2, j, p)^(1/6)*St*(ta – t(2, j, p))*(v(3, j, p) – v(1, j, p)) + 4*t(2, j, p)^2*(v(3, j, p) – 2*v(2, j, p) + v(1, j, p)) + 4*t(2, j, p)^2*v(2, j, p)*-2;
J11(1, 2) = -2*h*k*v(2, j, p)^(7/6)*St*(ta – t(2, j, p)) – t(2, j, p)^2*(2*v(3, j, p) – 2*v(1, j, p)) + 4*t(2, j, p)^2*v(2, j, p);
J11(end, end – 1) = 2*h*k*v(end – 1, j , p)^(7/6)*St*(ta – t(end – 1, j, p)) – t(end – 1, j , p)^2*(2*v(end – 2, j, p) – 2*v(end, j, p)) + 4*t(end – 1, j, p)^2*v(end – 1, j, p);
J11(end, end) = -2*h*k*(7/6)*v(end – 1, j, p)^(1/6)*St*(ta – t(end – 1, j , p))*(v(end, j, p) – v(end – 2, j, p)) + 4*t(end – 1, j, p)^2*(v(end, j, p) – 2*v(end – 1, j, p) + v(end – 2, j, p)) + 4*t(end – 1, j, p)^2*v(end – 1, j, p)*-2;
for i = 3 : 1 : m – 1
J11(i – 1, i – 2) = 2*h*k*v(i, j, p)^(7/6)*St*(ta – t(i, j, p)) – t(i)^2*(2*v(i – 1, j, p) – 2*v(i + 1, j, p)) + 4*t(i, j, p)^2*v(i, j, p);
J11(i – 1, i – 1) = -2*h*k*(7/6)*v(i, j, p)^(1/6)*St*(ta – t(i, j, p))*(v(i + 1, j, p) – v(i – 1, j, p)) + 4*t(i, j, p)^2*(v(i + 1, j, p) – 2*v(i, j, p) + v(i – 1, j, p)) + 4*t(i, j, p)^2*v(i, j, p)*-2;
J11(i – 1, i) = -2*h*k*v(i, j, p)^(7/6)*St*(ta – t(i, j, p)) – t(i, j, p)^2*(2*v(i + 1, j, p) – 2*v(i – 1, j, p)) + 4*t(i, j, p)^2*v(i, j, p);
end
J12 = zeros(m – 1, m);
J12(1, 1) = 2*h*k*v(2, j , p)^(7/6)*St*(v(3, j, p) – v(1, j, p)) – 2*t(2, j, p)*(v(3, j, p) – v(1, j, p))^2 + 8*t(2, j, p)*v(2, j, p)*(v(3, j, p) – 2*v(2, j, p) + v(1, j, p));
J12(end, end – 1) = 2*h*k*v(end – 1, j, p)^(7/6)*St*(v(end, j, p) – v(end – 2, j, p)) – 2*t(end – 1, j, p)*(v(end, j, p) – v(end – 2, j, p))^2 + 8*t(end – 1, j, p)*v(end – 1, j, p)*(v(end, j, p) – 2*v(end – 1, j, p) + v(end – 2, j, p));
for i = 3 : 1 : m – 1
J12(i – 1, i – 1) = 2*h*k*v(i, j, p)^(7/6)*St*(v(i + 1, j, p) – v(i – 1, j, p)) – 2*t(i, j, p)*(v(i + 1, j, p) – v(i – 1, j, p))^2 + 8*t(i, j, p)*v(i, j, p)*(v(i + 1, j, p) – 2*v(i, j, p) + v(i – 1, j, p));
end
J21 = zeros(m, m – 1);
for i = 2 : 1 : m
J21(i – 1, i – 1) = -2*h*St*(1/6)*v(i, j, p)^(-5/6)*(ta – t(i, j, p));
end
J22 = zeros(m, m);
J22(1, 1) = 2*h*St*v(2, j, p)^(1/6);
J22(1, 2) = 1;
J22(end, end – 2) = 1;
J22(end, end – 1) = -4;
J22(end, end) = 3 + 2*h*St*v(end, j, p)^(1/6);
for i = 3 : 1 : m
J22(i – 1, i – 2) = -1;
J22(i – 1, i – 1) = 2*h*St*v(i, j, p)^(1/6);
J22(i – 1, i) = 1;
end
J = [J11 J12 ; J21 J22];
%% Updating & Condition number
delta = -inv(J)*R;
J_i = inv(J);
for i = 1 : 1 : length(J)
aa(i) = sum(abs(J(i, :)));
bb(i) = sum(abs((J_i(i, :))));
end
aam = max(aa);
bbm = max(bb);
c = aam*bbm;
C(j, p) = 1/(aam*bbm);
%% Escaping condition & Condition number of jacobian0
if sqrt(sum(delta.^2)) <= tol
vf(:, p) = v(:, j, p);
tf(:, p) = t(:, j, p);
if p <= length(St_v) – 1
j = 1;
v(:, j, p + 1) = vf(:, p);
t(:, j, p + 1) = tf(:, p);
else
v(:, j, p + 1) = 0;
t(:, j, p + 1) = 0;
end
break
else
v(:, j + 1, p) = v(:, j, p);
v(2 : 1 : end – 1, j + 1, p) = v(2 : 1 : end – 1, j, p) + delta(1 : 1 : m – 1);
t(:, j + 1, p) = t(:, j, p);
t(2 : 1 : end, j + 1, p) = t(2 : 1 : end, j, p) + delta(m : 1 : end);
j = j + 1;
iter = iter + 1;
end
end
%% Plotting
plot(x, vf(:, p), ‘r’)
hold on
plot(x, tf(:, p), ‘b’)
hold on
end differential equations, newton-raphson method MATLAB Answers — New Questions