Category: Matlab
Category Archives: Matlab
i want matrix E == matrix B(i want correct answer)
I have tried to run this code many times , even changed the equation , values but i am not able to get the correct answer , here A given is before its transpose of glcm
% Convert to grayscale if necessary
%if size(img, 3) > 1
%img = rgb2gray(img);
%end
I = [0 0 1 1;
0 0 1 1;
0 2 2 2;
2 2 3 3];
A = [2 2 1 0;
0 2 0 0;
0 0 3 1;
0 0 0 1];
D = A’;
E = A + D;
E
%glcms = graycomatrix(I);
%glcms
m = 5;
%n = 4;
count = 0;
%count1 = 0;
n = 3;
o = 4;
for i = 1 : 4
for j = 1:4
for l = 1 : 4
for k = 1 : 4
x = I(i,1);
y = I(i,3-1);
a = I(i,3+1-1);
b = I(i,2+2-1);
if (I(i,j) && I(i,n-j)) == (I(i,n-k) && I(i, n+1-k))
%disp(‘Entered the control’);
% count = count+ 1;
%end
end
end
B(i,j) = count;
end
%B(i,j) = count;
%end
%end
end
%if (E == B)
% disp(‘Correct Answer’);
%else
% disp(‘Not Correct Answer’);
%endI have tried to run this code many times , even changed the equation , values but i am not able to get the correct answer , here A given is before its transpose of glcm
% Convert to grayscale if necessary
%if size(img, 3) > 1
%img = rgb2gray(img);
%end
I = [0 0 1 1;
0 0 1 1;
0 2 2 2;
2 2 3 3];
A = [2 2 1 0;
0 2 0 0;
0 0 3 1;
0 0 0 1];
D = A’;
E = A + D;
E
%glcms = graycomatrix(I);
%glcms
m = 5;
%n = 4;
count = 0;
%count1 = 0;
n = 3;
o = 4;
for i = 1 : 4
for j = 1:4
for l = 1 : 4
for k = 1 : 4
x = I(i,1);
y = I(i,3-1);
a = I(i,3+1-1);
b = I(i,2+2-1);
if (I(i,j) && I(i,n-j)) == (I(i,n-k) && I(i, n+1-k))
%disp(‘Entered the control’);
% count = count+ 1;
%end
end
end
B(i,j) = count;
end
%B(i,j) = count;
%end
%end
end
%if (E == B)
% disp(‘Correct Answer’);
%else
% disp(‘Not Correct Answer’);
%end I have tried to run this code many times , even changed the equation , values but i am not able to get the correct answer , here A given is before its transpose of glcm
% Convert to grayscale if necessary
%if size(img, 3) > 1
%img = rgb2gray(img);
%end
I = [0 0 1 1;
0 0 1 1;
0 2 2 2;
2 2 3 3];
A = [2 2 1 0;
0 2 0 0;
0 0 3 1;
0 0 0 1];
D = A’;
E = A + D;
E
%glcms = graycomatrix(I);
%glcms
m = 5;
%n = 4;
count = 0;
%count1 = 0;
n = 3;
o = 4;
for i = 1 : 4
for j = 1:4
for l = 1 : 4
for k = 1 : 4
x = I(i,1);
y = I(i,3-1);
a = I(i,3+1-1);
b = I(i,2+2-1);
if (I(i,j) && I(i,n-j)) == (I(i,n-k) && I(i, n+1-k))
%disp(‘Entered the control’);
% count = count+ 1;
%end
end
end
B(i,j) = count;
end
%B(i,j) = count;
%end
%end
end
%if (E == B)
% disp(‘Correct Answer’);
%else
% disp(‘Not Correct Answer’);
%end glcm MATLAB Answers — New Questions
Problem with Section 7.1 Task 5 of Simscape onramp
Although I correctly configured the model blocks and parameters accordingly to the directions, the training evaluation said my output did not match the answer. It prevents me from proceeding to the next task.
Here I attach two snapshots of the evaluation plot of Section 7.1 Task 5 where the problem occured.
The first picture shows comparison of the signals spanning to the end time (10sec). It looks to me exact agreement between the requirement and my signal, but it indicates slight discrepancies in the beginning of simulation (from 0 sec to around 1.0 sec) shown by the red dots. Due to these errors, the evaluator says my signal is incorrect.
The second picture zooms into the details of the red dots. It looks to me my signal has a slight delay to the answer.
I confirmed that my setup is consistent with the directions and the results unchanged even after several retries from the beginning of Section 7.1.
I am wondering if there is another parameter to be adjusted, which is not explicitly directed.
Can I have your help for solving this problem?Although I correctly configured the model blocks and parameters accordingly to the directions, the training evaluation said my output did not match the answer. It prevents me from proceeding to the next task.
Here I attach two snapshots of the evaluation plot of Section 7.1 Task 5 where the problem occured.
The first picture shows comparison of the signals spanning to the end time (10sec). It looks to me exact agreement between the requirement and my signal, but it indicates slight discrepancies in the beginning of simulation (from 0 sec to around 1.0 sec) shown by the red dots. Due to these errors, the evaluator says my signal is incorrect.
The second picture zooms into the details of the red dots. It looks to me my signal has a slight delay to the answer.
I confirmed that my setup is consistent with the directions and the results unchanged even after several retries from the beginning of Section 7.1.
I am wondering if there is another parameter to be adjusted, which is not explicitly directed.
Can I have your help for solving this problem? Although I correctly configured the model blocks and parameters accordingly to the directions, the training evaluation said my output did not match the answer. It prevents me from proceeding to the next task.
Here I attach two snapshots of the evaluation plot of Section 7.1 Task 5 where the problem occured.
The first picture shows comparison of the signals spanning to the end time (10sec). It looks to me exact agreement between the requirement and my signal, but it indicates slight discrepancies in the beginning of simulation (from 0 sec to around 1.0 sec) shown by the red dots. Due to these errors, the evaluator says my signal is incorrect.
The second picture zooms into the details of the red dots. It looks to me my signal has a slight delay to the answer.
I confirmed that my setup is consistent with the directions and the results unchanged even after several retries from the beginning of Section 7.1.
I am wondering if there is another parameter to be adjusted, which is not explicitly directed.
Can I have your help for solving this problem? simscape, onramp MATLAB Answers — New Questions
Error: The “A” matrix must be a numeric array with no Inf’s or NaN’s.
coder.extrinsic("testcon2dis")
I_2 = eye(2);
O_2 = zeros(2);
Ac_xy = [O_2 -(1/Lfilter)*I_2 O_2; (1/Cfilter)*I_2 O_2 -(1/Cfilter)*I_2; O_2 O_2 w* J];
Bc_xy = [(Vdc/Lfilter)*I_2 O_2 O_2]’ * Tr;
Cxy = [I_2 O_2 O_2; O_2 I_2 O_2];
[A, B] = testcon2dis(Ac_xy, Bc_xy, Cxy, Ts);
%%%%%%%%%%%%%%%% Below is the testcon2dis%%%%%%%%%%%%%%%%
function[A,B] = testcon2dis(Ac_xy,Bc_xy,Cxy,Ts)
ct_sys = ss(Ac_xy,Bc_xy,Cxy,[]);
dt_sys = c2d(ct_sys,Ts);
A = dt_sys.a;
B = dt_sys.b;
end
%%%%%%%%%%%%%%%%%%%%%%% Error message
Error:The "A" matrix must be a numeric array with no Inf’s or NaN’s.
Error in ss.ss.m (line 290)
throw(ME)
Error in testcon2dis.m (line 2)
Error in ‘MPC_3Phase_Inverter/MATLAB Function1’ (line 31)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Any suggestions I tried calling the A marix as a double or zeros to give it sizes but still got new errors doesn’t work(this block of code is used in a matlab function block simulink)
%%%%%%%%%%%%%%%%%%%%% Error message calling A as a double or zeros(6,6)%%%%%%%%%%%%%%%
Warning: Matrix is singular, close to singular or badly scaled. Results may be inaccurate. RCOND = NaN.
> In MPC_Matrices_l (line 69)
Error:Matrix must be positive definite.
Error in MPC_Matrices_l.m (line 70)
Hinv=chol(W_inv);
Error in ‘MPC_3Phase_Inverter/MATLAB Function1’ (line 43)
Any suggestions thank youcoder.extrinsic("testcon2dis")
I_2 = eye(2);
O_2 = zeros(2);
Ac_xy = [O_2 -(1/Lfilter)*I_2 O_2; (1/Cfilter)*I_2 O_2 -(1/Cfilter)*I_2; O_2 O_2 w* J];
Bc_xy = [(Vdc/Lfilter)*I_2 O_2 O_2]’ * Tr;
Cxy = [I_2 O_2 O_2; O_2 I_2 O_2];
[A, B] = testcon2dis(Ac_xy, Bc_xy, Cxy, Ts);
%%%%%%%%%%%%%%%% Below is the testcon2dis%%%%%%%%%%%%%%%%
function[A,B] = testcon2dis(Ac_xy,Bc_xy,Cxy,Ts)
ct_sys = ss(Ac_xy,Bc_xy,Cxy,[]);
dt_sys = c2d(ct_sys,Ts);
A = dt_sys.a;
B = dt_sys.b;
end
%%%%%%%%%%%%%%%%%%%%%%% Error message
Error:The "A" matrix must be a numeric array with no Inf’s or NaN’s.
Error in ss.ss.m (line 290)
throw(ME)
Error in testcon2dis.m (line 2)
Error in ‘MPC_3Phase_Inverter/MATLAB Function1’ (line 31)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Any suggestions I tried calling the A marix as a double or zeros to give it sizes but still got new errors doesn’t work(this block of code is used in a matlab function block simulink)
%%%%%%%%%%%%%%%%%%%%% Error message calling A as a double or zeros(6,6)%%%%%%%%%%%%%%%
Warning: Matrix is singular, close to singular or badly scaled. Results may be inaccurate. RCOND = NaN.
> In MPC_Matrices_l (line 69)
Error:Matrix must be positive definite.
Error in MPC_Matrices_l.m (line 70)
Hinv=chol(W_inv);
Error in ‘MPC_3Phase_Inverter/MATLAB Function1’ (line 43)
Any suggestions thank you coder.extrinsic("testcon2dis")
I_2 = eye(2);
O_2 = zeros(2);
Ac_xy = [O_2 -(1/Lfilter)*I_2 O_2; (1/Cfilter)*I_2 O_2 -(1/Cfilter)*I_2; O_2 O_2 w* J];
Bc_xy = [(Vdc/Lfilter)*I_2 O_2 O_2]’ * Tr;
Cxy = [I_2 O_2 O_2; O_2 I_2 O_2];
[A, B] = testcon2dis(Ac_xy, Bc_xy, Cxy, Ts);
%%%%%%%%%%%%%%%% Below is the testcon2dis%%%%%%%%%%%%%%%%
function[A,B] = testcon2dis(Ac_xy,Bc_xy,Cxy,Ts)
ct_sys = ss(Ac_xy,Bc_xy,Cxy,[]);
dt_sys = c2d(ct_sys,Ts);
A = dt_sys.a;
B = dt_sys.b;
end
%%%%%%%%%%%%%%%%%%%%%%% Error message
Error:The "A" matrix must be a numeric array with no Inf’s or NaN’s.
Error in ss.ss.m (line 290)
throw(ME)
Error in testcon2dis.m (line 2)
Error in ‘MPC_3Phase_Inverter/MATLAB Function1’ (line 31)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Any suggestions I tried calling the A marix as a double or zeros to give it sizes but still got new errors doesn’t work(this block of code is used in a matlab function block simulink)
%%%%%%%%%%%%%%%%%%%%% Error message calling A as a double or zeros(6,6)%%%%%%%%%%%%%%%
Warning: Matrix is singular, close to singular or badly scaled. Results may be inaccurate. RCOND = NaN.
> In MPC_Matrices_l (line 69)
Error:Matrix must be positive definite.
Error in MPC_Matrices_l.m (line 70)
Hinv=chol(W_inv);
Error in ‘MPC_3Phase_Inverter/MATLAB Function1’ (line 43)
Any suggestions thank you matlab, simulink, matlab function MATLAB Answers — New Questions
cross correlation between two time series data
I have two different normalised atmospheric time series data. Their temporal variation shows a mixed pattern of correlation and anti correlation at different time periods. Aparently there are some lags in these correlations. But when I carry out cross correlation, the plot shows max correlation at 1500 time lag which is the number of data points (plot xcorr-data.png). The wcoherence shows strong correlation between 0.5-2 days buts arrows dont show any definite pattern. Can anyone suggest me, how I can interprete the result in more meaningful way. Any suggestion about any alternative ways will be great.
Thanking you in advance.I have two different normalised atmospheric time series data. Their temporal variation shows a mixed pattern of correlation and anti correlation at different time periods. Aparently there are some lags in these correlations. But when I carry out cross correlation, the plot shows max correlation at 1500 time lag which is the number of data points (plot xcorr-data.png). The wcoherence shows strong correlation between 0.5-2 days buts arrows dont show any definite pattern. Can anyone suggest me, how I can interprete the result in more meaningful way. Any suggestion about any alternative ways will be great.
Thanking you in advance. I have two different normalised atmospheric time series data. Their temporal variation shows a mixed pattern of correlation and anti correlation at different time periods. Aparently there are some lags in these correlations. But when I carry out cross correlation, the plot shows max correlation at 1500 time lag which is the number of data points (plot xcorr-data.png). The wcoherence shows strong correlation between 0.5-2 days buts arrows dont show any definite pattern. Can anyone suggest me, how I can interprete the result in more meaningful way. Any suggestion about any alternative ways will be great.
Thanking you in advance. xcorr, wcoherence, time series, cross correlation, atmospheric data MATLAB Answers — New Questions
Help creating RMS Window
Hello,
I am supposed to import data from 10 Excel files then filter, rectify and create a RMS window for the data, but I am stuck at the last step. I am not very good at MATLAB, but I believe that I am on the right track so far. Any help would be extremely appreciated. I included what I have so far; the commented section at the end is what I am supposed to work off of. I tried a similar loop like the ones prior but did not succeed. Thank you in advance!
clear all; clc
fileName ={‘Gait_normal01.csv’;’Gait_normal02.csv’;’Gait_rightlimp01.csv’;’Gait_rightlimp02.csv’;’MVC_Extension01.csv’;…
‘MVC_Extension02.csv’;’MVC_Flexion01.csv’;’MVC_Flexion02.csv’;’SitStand01.csv’;’Squat01.csv’}
% read in Data
for i=1:10
EMGraw_data{i,1} =dlmread(fileName{i,1},’,’,5,2);
end
%filter
for i=1:10
i
for j=1:4
j
[B,A]=butter(2,[10/500 350/500],’bandpass’);
EMGf{i,1}(:,j) = filtfilt(B,A,EMGraw_data{i,1}(:,j));
end
end
% rectify
for i=1:10
i
for j=1:4
j
EMGrec{i,1}(:,j) = abs(EMGf{i,1}(:,j));
end
end
% %Root mean square
% winlength = 299; %input 1-desired window length
% EMGrms=zeros(length(EMGrec)-winlength,1);
% for i=1:length(EMGrec)-winlength;
% win=EMGrec(i:i+winlength);
% EMGrms(i+(winlength+1)/2)=sqrt(sum(win.^2)/winlength);
% endHello,
I am supposed to import data from 10 Excel files then filter, rectify and create a RMS window for the data, but I am stuck at the last step. I am not very good at MATLAB, but I believe that I am on the right track so far. Any help would be extremely appreciated. I included what I have so far; the commented section at the end is what I am supposed to work off of. I tried a similar loop like the ones prior but did not succeed. Thank you in advance!
clear all; clc
fileName ={‘Gait_normal01.csv’;’Gait_normal02.csv’;’Gait_rightlimp01.csv’;’Gait_rightlimp02.csv’;’MVC_Extension01.csv’;…
‘MVC_Extension02.csv’;’MVC_Flexion01.csv’;’MVC_Flexion02.csv’;’SitStand01.csv’;’Squat01.csv’}
% read in Data
for i=1:10
EMGraw_data{i,1} =dlmread(fileName{i,1},’,’,5,2);
end
%filter
for i=1:10
i
for j=1:4
j
[B,A]=butter(2,[10/500 350/500],’bandpass’);
EMGf{i,1}(:,j) = filtfilt(B,A,EMGraw_data{i,1}(:,j));
end
end
% rectify
for i=1:10
i
for j=1:4
j
EMGrec{i,1}(:,j) = abs(EMGf{i,1}(:,j));
end
end
% %Root mean square
% winlength = 299; %input 1-desired window length
% EMGrms=zeros(length(EMGrec)-winlength,1);
% for i=1:length(EMGrec)-winlength;
% win=EMGrec(i:i+winlength);
% EMGrms(i+(winlength+1)/2)=sqrt(sum(win.^2)/winlength);
% end Hello,
I am supposed to import data from 10 Excel files then filter, rectify and create a RMS window for the data, but I am stuck at the last step. I am not very good at MATLAB, but I believe that I am on the right track so far. Any help would be extremely appreciated. I included what I have so far; the commented section at the end is what I am supposed to work off of. I tried a similar loop like the ones prior but did not succeed. Thank you in advance!
clear all; clc
fileName ={‘Gait_normal01.csv’;’Gait_normal02.csv’;’Gait_rightlimp01.csv’;’Gait_rightlimp02.csv’;’MVC_Extension01.csv’;…
‘MVC_Extension02.csv’;’MVC_Flexion01.csv’;’MVC_Flexion02.csv’;’SitStand01.csv’;’Squat01.csv’}
% read in Data
for i=1:10
EMGraw_data{i,1} =dlmread(fileName{i,1},’,’,5,2);
end
%filter
for i=1:10
i
for j=1:4
j
[B,A]=butter(2,[10/500 350/500],’bandpass’);
EMGf{i,1}(:,j) = filtfilt(B,A,EMGraw_data{i,1}(:,j));
end
end
% rectify
for i=1:10
i
for j=1:4
j
EMGrec{i,1}(:,j) = abs(EMGf{i,1}(:,j));
end
end
% %Root mean square
% winlength = 299; %input 1-desired window length
% EMGrms=zeros(length(EMGrec)-winlength,1);
% for i=1:length(EMGrec)-winlength;
% win=EMGrec(i:i+winlength);
% EMGrms(i+(winlength+1)/2)=sqrt(sum(win.^2)/winlength);
% end rms, for, loop, for loop, win, winlength, root, mean, square MATLAB Answers — New Questions
Stateful predict block failed to evaluate mask initial commands
this is the error showing up as soon as I run my model.
The recurrent neural network that I’ve used in this stateful predict block works completely fine on its own but when I use its MAT file in this block.. it ain’t working.
what should I dothis is the error showing up as soon as I run my model.
The recurrent neural network that I’ve used in this stateful predict block works completely fine on its own but when I use its MAT file in this block.. it ain’t working.
what should I do this is the error showing up as soon as I run my model.
The recurrent neural network that I’ve used in this stateful predict block works completely fine on its own but when I use its MAT file in this block.. it ain’t working.
what should I do stateful predict, simulink, simulation, deep learning, neural network MATLAB Answers — New Questions
How to generate code for self attention layer for deployment
How to generate code for self attention layer for deploymentHow to generate code for self attention layer for deployment How to generate code for self attention layer for deployment code generation MATLAB Answers — New Questions
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
error using sym/solve and maplemex
Hi! Whenever I try and solve a systems of equations such as the simple one above, an error like this one shows up, where it says that there was a problem using sym/solve and maplemex. This is a problem I’ve had since December and I have no idea how to fix it. I have the symbolic math toolbox installed, so that does not seem to be the issue. Any help would be greatly appreciated!Hi! Whenever I try and solve a systems of equations such as the simple one above, an error like this one shows up, where it says that there was a problem using sym/solve and maplemex. This is a problem I’ve had since December and I have no idea how to fix it. I have the symbolic math toolbox installed, so that does not seem to be the issue. Any help would be greatly appreciated! Hi! Whenever I try and solve a systems of equations such as the simple one above, an error like this one shows up, where it says that there was a problem using sym/solve and maplemex. This is a problem I’ve had since December and I have no idea how to fix it. I have the symbolic math toolbox installed, so that does not seem to be the issue. Any help would be greatly appreciated! #syms #error #toolbox #solve #systems #of #equations #maplemex #symbolic, symbolic, syms, solve MATLAB Answers — New Questions
How to mark step changes of various sizes
This is my code for part of a project im doing:
time_run2 = transpose((1:numel(second_run_data(:, 1))’) / 60);
figure(3);
plot(time_run2, second_run_data(:, 1))
whish is used to plot a section of data from a larger set. This data follows a step pattern however each step is brocken up into a bunch of little steps and then a section of flat followed by another step and more flat. some of the steps are slightly curved so i cant simply use the ischange() function to find them becasue it wil find multiple changes in each step. My task is too mark each step with a circle or an x whatever but im only allowed to have 1 mark per step. Is there any way of doing this becasue I have no ideaThis is my code for part of a project im doing:
time_run2 = transpose((1:numel(second_run_data(:, 1))’) / 60);
figure(3);
plot(time_run2, second_run_data(:, 1))
whish is used to plot a section of data from a larger set. This data follows a step pattern however each step is brocken up into a bunch of little steps and then a section of flat followed by another step and more flat. some of the steps are slightly curved so i cant simply use the ischange() function to find them becasue it wil find multiple changes in each step. My task is too mark each step with a circle or an x whatever but im only allowed to have 1 mark per step. Is there any way of doing this becasue I have no idea This is my code for part of a project im doing:
time_run2 = transpose((1:numel(second_run_data(:, 1))’) / 60);
figure(3);
plot(time_run2, second_run_data(:, 1))
whish is used to plot a section of data from a larger set. This data follows a step pattern however each step is brocken up into a bunch of little steps and then a section of flat followed by another step and more flat. some of the steps are slightly curved so i cant simply use the ischange() function to find them becasue it wil find multiple changes in each step. My task is too mark each step with a circle or an x whatever but im only allowed to have 1 mark per step. Is there any way of doing this becasue I have no idea step, data, plot MATLAB Answers — New Questions
lost in my raw unstructured data files, don’t know where is my foot stand , and keep rounding in cycle
envision a case where a user with no code exp , or data analysis exp , found him self in the peak of AI trend , he followed his dream, and passion to building his project system software , and for 15 months, he just into AI playground , from researcher , experiments , codes , concepts , ideas brain storming , think tanks , every day and night, just running randomly , without plan , or road map , now after this much time he know he reach some thing big , but when he stop to embody and embark it, he found him self stuck in multi raw unstructured data files , unorganized , don;t know where his foot stand , now he keep round in cycle without knowing how to move forword , now he have this idea as all this unstructured data saved in google drive , he want to use google colab jypter notebook , + google gemnini 1.5 API , + langchain , lang graph , llama index , to clean and structured this data into knowledge graph database , where he can analysis every thing connect it , find the hidden pattern, realationship between the data , or any other ideas , so he can find where his foot stand , and how he can move forword, im this user 🙁 , and really need helpenvision a case where a user with no code exp , or data analysis exp , found him self in the peak of AI trend , he followed his dream, and passion to building his project system software , and for 15 months, he just into AI playground , from researcher , experiments , codes , concepts , ideas brain storming , think tanks , every day and night, just running randomly , without plan , or road map , now after this much time he know he reach some thing big , but when he stop to embody and embark it, he found him self stuck in multi raw unstructured data files , unorganized , don;t know where his foot stand , now he keep round in cycle without knowing how to move forword , now he have this idea as all this unstructured data saved in google drive , he want to use google colab jypter notebook , + google gemnini 1.5 API , + langchain , lang graph , llama index , to clean and structured this data into knowledge graph database , where he can analysis every thing connect it , find the hidden pattern, realationship between the data , or any other ideas , so he can find where his foot stand , and how he can move forword, im this user 🙁 , and really need help envision a case where a user with no code exp , or data analysis exp , found him self in the peak of AI trend , he followed his dream, and passion to building his project system software , and for 15 months, he just into AI playground , from researcher , experiments , codes , concepts , ideas brain storming , think tanks , every day and night, just running randomly , without plan , or road map , now after this much time he know he reach some thing big , but when he stop to embody and embark it, he found him self stuck in multi raw unstructured data files , unorganized , don;t know where his foot stand , now he keep round in cycle without knowing how to move forword , now he have this idea as all this unstructured data saved in google drive , he want to use google colab jypter notebook , + google gemnini 1.5 API , + langchain , lang graph , llama index , to clean and structured this data into knowledge graph database , where he can analysis every thing connect it , find the hidden pattern, realationship between the data , or any other ideas , so he can find where his foot stand , and how he can move forword, im this user 🙁 , and really need help unstructured data, database, knowledge graph MATLAB Answers — New Questions
Calculation of centroide of an object in image
Hello everyone
i am currently working on a school project of medical image. I have an x-ray image of a head in lateral position with a brain tumor. i have to calculate the centroide of this tumor. the tumor has an higher contrast however so does other objects in the image. how can i do it?
thank you everyone in advance.
i don’t have much experience with matlab so i am a bit in a struggleHello everyone
i am currently working on a school project of medical image. I have an x-ray image of a head in lateral position with a brain tumor. i have to calculate the centroide of this tumor. the tumor has an higher contrast however so does other objects in the image. how can i do it?
thank you everyone in advance.
i don’t have much experience with matlab so i am a bit in a struggle Hello everyone
i am currently working on a school project of medical image. I have an x-ray image of a head in lateral position with a brain tumor. i have to calculate the centroide of this tumor. the tumor has an higher contrast however so does other objects in the image. how can i do it?
thank you everyone in advance.
i don’t have much experience with matlab so i am a bit in a struggle medical image, image analysis MATLAB Answers — New Questions
How to fuse multiple images while retaining the color of the original images and changing some image colors?
Hello, all,
I am trying to fuse multiple images together to form a final image with the original colors in the original images, however, when I fuse them together, they are mixing and blending the colors to be a completely different color from what I am looking for.
Can I get some insight on solving this please?
Final Image that needs help fixing the colors:
Full map image stored as 600x600x3 uint8:
Fastest path image (the white pixels here need to be green to indicate the fastest path in the final image) stored as 600x600x3 uint8:
Trajectory recorded (the black and magenta pixels) not the exact example from the above images but they are in the same structure:
Stored in workspace as: matImgRGB in a 600×600 uint8 and this is the outcome when using infuse() with the Fastest path image above:
The final image should be the green fastest path image with the colors of the tracjectory recorded overlaid with the full map image with it’s blue color. Appreciate the help!
Many thanks,
MichaelHello, all,
I am trying to fuse multiple images together to form a final image with the original colors in the original images, however, when I fuse them together, they are mixing and blending the colors to be a completely different color from what I am looking for.
Can I get some insight on solving this please?
Final Image that needs help fixing the colors:
Full map image stored as 600x600x3 uint8:
Fastest path image (the white pixels here need to be green to indicate the fastest path in the final image) stored as 600x600x3 uint8:
Trajectory recorded (the black and magenta pixels) not the exact example from the above images but they are in the same structure:
Stored in workspace as: matImgRGB in a 600×600 uint8 and this is the outcome when using infuse() with the Fastest path image above:
The final image should be the green fastest path image with the colors of the tracjectory recorded overlaid with the full map image with it’s blue color. Appreciate the help!
Many thanks,
Michael Hello, all,
I am trying to fuse multiple images together to form a final image with the original colors in the original images, however, when I fuse them together, they are mixing and blending the colors to be a completely different color from what I am looking for.
Can I get some insight on solving this please?
Final Image that needs help fixing the colors:
Full map image stored as 600x600x3 uint8:
Fastest path image (the white pixels here need to be green to indicate the fastest path in the final image) stored as 600x600x3 uint8:
Trajectory recorded (the black and magenta pixels) not the exact example from the above images but they are in the same structure:
Stored in workspace as: matImgRGB in a 600×600 uint8 and this is the outcome when using infuse() with the Fastest path image above:
The final image should be the green fastest path image with the colors of the tracjectory recorded overlaid with the full map image with it’s blue color. Appreciate the help!
Many thanks,
Michael image processing MATLAB Answers — New Questions
How to increase Buffer size parameter in MATLAB .mdl simulation file?
I am having error like "The specified buffer was too small.During simulation, the buffer size was temporarily increased to 2048. In order to use Real-Time Workshop, you need to update the buffer size parameter"I am having error like "The specified buffer was too small.During simulation, the buffer size was temporarily increased to 2048. In order to use Real-Time Workshop, you need to update the buffer size parameter" I am having error like "The specified buffer was too small.During simulation, the buffer size was temporarily increased to 2048. In order to use Real-Time Workshop, you need to update the buffer size parameter" figure, circuit MATLAB Answers — New Questions
Not enough input arguments ode45
Trying to find theta, angular velocity, and angular acceleration.
clear all
close all
clc
syms theta1(t)
%physical constants
Dd=0.1; Dw=0.75; Dh=0.5; Dm=15; ax=0.1; ay=0.5; b=0.25;
a = sqrt(ax^2 + ay^2); Jd = (1/3)*Dm*Dh^2;
Dweight=Dm*9.81;
thetaD0 = 0;
dthetaDdt0=0;
PHI = atand(ax/ay);
Lo = sqrt(a^2 + b^2 -2*a*b*cosd(PHI));
xmax = sqrt(a^2 + b^2 -2*a*b*cosd(PHI+90)) – Lo;
Ltotal = Lo + xmax;
k=500;
numphi = a*sind(PHI + thetaD0);
denphi = sqrt(a^2 + b^2 -2*a*b*cosd(PHI+thetaD0));
phi = asind(numphi/denphi);
d = 0.01;
%motor
NoLoadSpeed = 1057.672; %rad/s
Kt = 0.01386; %Vsec/rad
NoLoadI = 0.036; %A
Cvf = Kt*NoLoadI/NoLoadSpeed;
Jm =4.2E-7; %kg m^2
w = 1.5708/4 ; %rad/s
x = sqrt(a^2 + b^2 -2*a*b*cosd(PHI+theta1)) – Lo;
i=0.036;
[t, sol] = calcsum(thetaD0, dthetaDdt0, Cvf, b, phi, d, Dh, Dm, k, i, Kt, x, Jd, Jm)
%sol = [0;0;0]
%derivs = myODE(t, sol, Cvf, b, phi, d, Dh, Dm, k, i, Kt, x, thetaD0, Jd, Jm)
%%
function [t, sol] = calcsum(thetaD0, dthetaDdt0, Cvf, b, phi, d, Dh, Dm, k, i, Kt, x, Jd, Jm)
sol0 = [thetaD0; dthetaDdt0; i];
tspan = [0, 5];
[t, sol] = ode45(@(t, sol) myODE(sol, Cvf, b, phi, d, Dh, Dm, k, i, Kt, x, thetaD0, Jd, Jm), tspan, sol0);
end
%%
function derivs = myODE(t, sol, Cvf, b, phi, d, Dh, Dm, k, i, Kt, x, thetaD, Jd, Jm) %calc
dthetaDdt = sol(1);
g = (Kt*i – Cvf*dthetaDdt*(2*pi*b*sind(phi)/d));
num = 2*((2*pi/d)*((Kt*i – Cvf*dthetaDdt*(2*pi*b*sind(phi)/d))) + k*x)*b*sind(phi) – 0.5*Dm*9.81*Dh*sin(thetaD);
h = pi*b*sin(phi)/d;
den = Jd + 8*Jm*h^2;
ddthetaDdt = num/den;
derivs = [dthetaDdt; ddthetaDdt];
endTrying to find theta, angular velocity, and angular acceleration.
clear all
close all
clc
syms theta1(t)
%physical constants
Dd=0.1; Dw=0.75; Dh=0.5; Dm=15; ax=0.1; ay=0.5; b=0.25;
a = sqrt(ax^2 + ay^2); Jd = (1/3)*Dm*Dh^2;
Dweight=Dm*9.81;
thetaD0 = 0;
dthetaDdt0=0;
PHI = atand(ax/ay);
Lo = sqrt(a^2 + b^2 -2*a*b*cosd(PHI));
xmax = sqrt(a^2 + b^2 -2*a*b*cosd(PHI+90)) – Lo;
Ltotal = Lo + xmax;
k=500;
numphi = a*sind(PHI + thetaD0);
denphi = sqrt(a^2 + b^2 -2*a*b*cosd(PHI+thetaD0));
phi = asind(numphi/denphi);
d = 0.01;
%motor
NoLoadSpeed = 1057.672; %rad/s
Kt = 0.01386; %Vsec/rad
NoLoadI = 0.036; %A
Cvf = Kt*NoLoadI/NoLoadSpeed;
Jm =4.2E-7; %kg m^2
w = 1.5708/4 ; %rad/s
x = sqrt(a^2 + b^2 -2*a*b*cosd(PHI+theta1)) – Lo;
i=0.036;
[t, sol] = calcsum(thetaD0, dthetaDdt0, Cvf, b, phi, d, Dh, Dm, k, i, Kt, x, Jd, Jm)
%sol = [0;0;0]
%derivs = myODE(t, sol, Cvf, b, phi, d, Dh, Dm, k, i, Kt, x, thetaD0, Jd, Jm)
%%
function [t, sol] = calcsum(thetaD0, dthetaDdt0, Cvf, b, phi, d, Dh, Dm, k, i, Kt, x, Jd, Jm)
sol0 = [thetaD0; dthetaDdt0; i];
tspan = [0, 5];
[t, sol] = ode45(@(t, sol) myODE(sol, Cvf, b, phi, d, Dh, Dm, k, i, Kt, x, thetaD0, Jd, Jm), tspan, sol0);
end
%%
function derivs = myODE(t, sol, Cvf, b, phi, d, Dh, Dm, k, i, Kt, x, thetaD, Jd, Jm) %calc
dthetaDdt = sol(1);
g = (Kt*i – Cvf*dthetaDdt*(2*pi*b*sind(phi)/d));
num = 2*((2*pi/d)*((Kt*i – Cvf*dthetaDdt*(2*pi*b*sind(phi)/d))) + k*x)*b*sind(phi) – 0.5*Dm*9.81*Dh*sin(thetaD);
h = pi*b*sin(phi)/d;
den = Jd + 8*Jm*h^2;
ddthetaDdt = num/den;
derivs = [dthetaDdt; ddthetaDdt];
end Trying to find theta, angular velocity, and angular acceleration.
clear all
close all
clc
syms theta1(t)
%physical constants
Dd=0.1; Dw=0.75; Dh=0.5; Dm=15; ax=0.1; ay=0.5; b=0.25;
a = sqrt(ax^2 + ay^2); Jd = (1/3)*Dm*Dh^2;
Dweight=Dm*9.81;
thetaD0 = 0;
dthetaDdt0=0;
PHI = atand(ax/ay);
Lo = sqrt(a^2 + b^2 -2*a*b*cosd(PHI));
xmax = sqrt(a^2 + b^2 -2*a*b*cosd(PHI+90)) – Lo;
Ltotal = Lo + xmax;
k=500;
numphi = a*sind(PHI + thetaD0);
denphi = sqrt(a^2 + b^2 -2*a*b*cosd(PHI+thetaD0));
phi = asind(numphi/denphi);
d = 0.01;
%motor
NoLoadSpeed = 1057.672; %rad/s
Kt = 0.01386; %Vsec/rad
NoLoadI = 0.036; %A
Cvf = Kt*NoLoadI/NoLoadSpeed;
Jm =4.2E-7; %kg m^2
w = 1.5708/4 ; %rad/s
x = sqrt(a^2 + b^2 -2*a*b*cosd(PHI+theta1)) – Lo;
i=0.036;
[t, sol] = calcsum(thetaD0, dthetaDdt0, Cvf, b, phi, d, Dh, Dm, k, i, Kt, x, Jd, Jm)
%sol = [0;0;0]
%derivs = myODE(t, sol, Cvf, b, phi, d, Dh, Dm, k, i, Kt, x, thetaD0, Jd, Jm)
%%
function [t, sol] = calcsum(thetaD0, dthetaDdt0, Cvf, b, phi, d, Dh, Dm, k, i, Kt, x, Jd, Jm)
sol0 = [thetaD0; dthetaDdt0; i];
tspan = [0, 5];
[t, sol] = ode45(@(t, sol) myODE(sol, Cvf, b, phi, d, Dh, Dm, k, i, Kt, x, thetaD0, Jd, Jm), tspan, sol0);
end
%%
function derivs = myODE(t, sol, Cvf, b, phi, d, Dh, Dm, k, i, Kt, x, thetaD, Jd, Jm) %calc
dthetaDdt = sol(1);
g = (Kt*i – Cvf*dthetaDdt*(2*pi*b*sind(phi)/d));
num = 2*((2*pi/d)*((Kt*i – Cvf*dthetaDdt*(2*pi*b*sind(phi)/d))) + k*x)*b*sind(phi) – 0.5*Dm*9.81*Dh*sin(thetaD);
h = pi*b*sin(phi)/d;
den = Jd + 8*Jm*h^2;
ddthetaDdt = num/den;
derivs = [dthetaDdt; ddthetaDdt];
end ode45 MATLAB Answers — New Questions
Please need help applying EEG seizure detection code on edf files instead of txt
Hello, I’ve been trying to test different available codes to understand how to design an EEG seizure detection algorithm.
I would like to know how to change the code here: https://www.mathworks.com/help/wavelet/ug/time-frequency-convolutional-network-for-eeg-data-classification.html so I can use it fir the CHB MIT data (edf files) instead of txt files from the Bonn data set if possible. The section below in particular:
ds = tabularTextDatastore(dataDir,"IncludeSubfolders",true,"FileExtensions",".txt");
extraTXT = contains(tds.Files,"__MACOSX");
tds.Files(extraTXT) = [];
labels = filenames2labels(tds.Files,"ExtractBetween",[1 1]);
ii = 1;
eegData = cell(numel(labels),1);
while hasdata(tds)
tsTable = read(tds);
ts = tsTable.Var1;
eegData{ii} = reshape(ts,1,[]);
ii = ii+1;
endHello, I’ve been trying to test different available codes to understand how to design an EEG seizure detection algorithm.
I would like to know how to change the code here: https://www.mathworks.com/help/wavelet/ug/time-frequency-convolutional-network-for-eeg-data-classification.html so I can use it fir the CHB MIT data (edf files) instead of txt files from the Bonn data set if possible. The section below in particular:
ds = tabularTextDatastore(dataDir,"IncludeSubfolders",true,"FileExtensions",".txt");
extraTXT = contains(tds.Files,"__MACOSX");
tds.Files(extraTXT) = [];
labels = filenames2labels(tds.Files,"ExtractBetween",[1 1]);
ii = 1;
eegData = cell(numel(labels),1);
while hasdata(tds)
tsTable = read(tds);
ts = tsTable.Var1;
eegData{ii} = reshape(ts,1,[]);
ii = ii+1;
end Hello, I’ve been trying to test different available codes to understand how to design an EEG seizure detection algorithm.
I would like to know how to change the code here: https://www.mathworks.com/help/wavelet/ug/time-frequency-convolutional-network-for-eeg-data-classification.html so I can use it fir the CHB MIT data (edf files) instead of txt files from the Bonn data set if possible. The section below in particular:
ds = tabularTextDatastore(dataDir,"IncludeSubfolders",true,"FileExtensions",".txt");
extraTXT = contains(tds.Files,"__MACOSX");
tds.Files(extraTXT) = [];
labels = filenames2labels(tds.Files,"ExtractBetween",[1 1]);
ii = 1;
eegData = cell(numel(labels),1);
while hasdata(tds)
tsTable = read(tds);
ts = tsTable.Var1;
eegData{ii} = reshape(ts,1,[]);
ii = ii+1;
end deep learning, eeg, eeg detection, edf files, neural network, matlab MATLAB Answers — New Questions
Dutation type doesn’t preserve seconds accuracy
I’ve been working on some time conversion routines and ran into an issue I didn’t expect. Basically, the duration type doesn’t preserve seconds accuracy. This is caused by the fact that the duration type stores the duration as a single variable called millis (i.e., milliseconds) instead of separate d,h,m,ms or h,m,ms. Why this would be coded in that way does not make any sense to me. E.g., because of this you can get issues such as this:
format longg
dt1 = datetime(2000,1,1,0,0,1.2345678912345);
disp(dt1.Second)
dt2 = datetime(2000,1,1);
[~,~,s] = hms(dt1-dt2)
All well and good, the seconds is represented to full precision in the datetime variable, and the precision of the seconds difference is preserved in the subtraction because the datetimes are close to one another. But if they are not:
dt3 = datetime(1900,1,1);
[~,~,s] = hms(dt1-dt3)
You suddenly get a different result because the difference got mashed into a single millis variable and that large days value caused precision loss in the seconds. And simply constructing a duration type from scratch doesn’t solve anything either because of the internal duration storage as a single millis variable. E.g.,
[~,~,s] = hms(duration(0,0,1.2345678912345))
[~,~,s] = hms(duration(1e6,0,1.2345678912345))
I know, I know … "You shouldn’t depend on this level of accuracy in floating point calculations" etc. But that’s not the point. The point is as long as you have a datetime type that stores individual d,h,m,s properties, then IMHO the duration type should follow suit and preserve accuracy in calculations as much as possible to match it. If the duration type stored the data internally as d,h,m,s (or ms or ns) and took care when doing the internal arithmetic, the above differences could be avoided. In fact, since the millis is a private internal variable, I think TMW could still do this in future versions of MATLAB without breaking any user code since there shouldn’t be any user code that accesses the millis property. Any chance of this, TMW?
And speaking of the datetime type, IMHO the internal variable should have been millis or nanos instead of seconds. That way standard Terrestrial Time epochs such as J2000 (January 1, 2000, 11:58:55.816 UTC) etc. could be represented exactly internally. But since datetime y,m,d,h,m,s properties are public this unfortunately can’t be changed now.I’ve been working on some time conversion routines and ran into an issue I didn’t expect. Basically, the duration type doesn’t preserve seconds accuracy. This is caused by the fact that the duration type stores the duration as a single variable called millis (i.e., milliseconds) instead of separate d,h,m,ms or h,m,ms. Why this would be coded in that way does not make any sense to me. E.g., because of this you can get issues such as this:
format longg
dt1 = datetime(2000,1,1,0,0,1.2345678912345);
disp(dt1.Second)
dt2 = datetime(2000,1,1);
[~,~,s] = hms(dt1-dt2)
All well and good, the seconds is represented to full precision in the datetime variable, and the precision of the seconds difference is preserved in the subtraction because the datetimes are close to one another. But if they are not:
dt3 = datetime(1900,1,1);
[~,~,s] = hms(dt1-dt3)
You suddenly get a different result because the difference got mashed into a single millis variable and that large days value caused precision loss in the seconds. And simply constructing a duration type from scratch doesn’t solve anything either because of the internal duration storage as a single millis variable. E.g.,
[~,~,s] = hms(duration(0,0,1.2345678912345))
[~,~,s] = hms(duration(1e6,0,1.2345678912345))
I know, I know … "You shouldn’t depend on this level of accuracy in floating point calculations" etc. But that’s not the point. The point is as long as you have a datetime type that stores individual d,h,m,s properties, then IMHO the duration type should follow suit and preserve accuracy in calculations as much as possible to match it. If the duration type stored the data internally as d,h,m,s (or ms or ns) and took care when doing the internal arithmetic, the above differences could be avoided. In fact, since the millis is a private internal variable, I think TMW could still do this in future versions of MATLAB without breaking any user code since there shouldn’t be any user code that accesses the millis property. Any chance of this, TMW?
And speaking of the datetime type, IMHO the internal variable should have been millis or nanos instead of seconds. That way standard Terrestrial Time epochs such as J2000 (January 1, 2000, 11:58:55.816 UTC) etc. could be represented exactly internally. But since datetime y,m,d,h,m,s properties are public this unfortunately can’t be changed now. I’ve been working on some time conversion routines and ran into an issue I didn’t expect. Basically, the duration type doesn’t preserve seconds accuracy. This is caused by the fact that the duration type stores the duration as a single variable called millis (i.e., milliseconds) instead of separate d,h,m,ms or h,m,ms. Why this would be coded in that way does not make any sense to me. E.g., because of this you can get issues such as this:
format longg
dt1 = datetime(2000,1,1,0,0,1.2345678912345);
disp(dt1.Second)
dt2 = datetime(2000,1,1);
[~,~,s] = hms(dt1-dt2)
All well and good, the seconds is represented to full precision in the datetime variable, and the precision of the seconds difference is preserved in the subtraction because the datetimes are close to one another. But if they are not:
dt3 = datetime(1900,1,1);
[~,~,s] = hms(dt1-dt3)
You suddenly get a different result because the difference got mashed into a single millis variable and that large days value caused precision loss in the seconds. And simply constructing a duration type from scratch doesn’t solve anything either because of the internal duration storage as a single millis variable. E.g.,
[~,~,s] = hms(duration(0,0,1.2345678912345))
[~,~,s] = hms(duration(1e6,0,1.2345678912345))
I know, I know … "You shouldn’t depend on this level of accuracy in floating point calculations" etc. But that’s not the point. The point is as long as you have a datetime type that stores individual d,h,m,s properties, then IMHO the duration type should follow suit and preserve accuracy in calculations as much as possible to match it. If the duration type stored the data internally as d,h,m,s (or ms or ns) and took care when doing the internal arithmetic, the above differences could be avoided. In fact, since the millis is a private internal variable, I think TMW could still do this in future versions of MATLAB without breaking any user code since there shouldn’t be any user code that accesses the millis property. Any chance of this, TMW?
And speaking of the datetime type, IMHO the internal variable should have been millis or nanos instead of seconds. That way standard Terrestrial Time epochs such as J2000 (January 1, 2000, 11:58:55.816 UTC) etc. could be represented exactly internally. But since datetime y,m,d,h,m,s properties are public this unfortunately can’t be changed now. duration, datetime, seconds MATLAB Answers — New Questions
Error due to not enough input arguments in ‘ C_c = cold.m_dot_c*Cold.c_p_c; %W/K ‘
function [ q, e, UA, T_c_o ] = find_UA( hot, cold )
%Determine minimum heat capacity rate
C_c = cold.m_dot_c * cold.c_p_c; %W/K
C_h = hot.m_dot_h * hot.c_p_h; %W/K
C_min = min(C_c, C_h);
C_max = max(C_c, C_h);function [ q, e, UA, T_c_o ] = find_UA( hot, cold )
%Determine minimum heat capacity rate
C_c = cold.m_dot_c * cold.c_p_c; %W/K
C_h = hot.m_dot_h * hot.c_p_h; %W/K
C_min = min(C_c, C_h);
C_max = max(C_c, C_h); function [ q, e, UA, T_c_o ] = find_UA( hot, cold )
%Determine minimum heat capacity rate
C_c = cold.m_dot_c * cold.c_p_c; %W/K
C_h = hot.m_dot_h * hot.c_p_h; %W/K
C_min = min(C_c, C_h);
C_max = max(C_c, C_h); input arguments MATLAB Answers — New Questions
How to modify the response of a antenna element with respect to different opearting frequency in the specified frequency range?
Take the phased.IsotropoicAntennaElement as an example. ‘FrequencyRange’ is set between 1GHz and 2GHz. Under the default cofiguration, the response of the antenna in the specified ‘FrequencyRange’ is a constant of one. I wonder how can I change the respose with respect to different frequency.
antenna = phased.IsotropicAntennaElement(‘FrequencyRange’,[1e9 2e9]);
fc = 1e9*(1:0.1:2);
resp = antenna(fc,[0;0])Take the phased.IsotropoicAntennaElement as an example. ‘FrequencyRange’ is set between 1GHz and 2GHz. Under the default cofiguration, the response of the antenna in the specified ‘FrequencyRange’ is a constant of one. I wonder how can I change the respose with respect to different frequency.
antenna = phased.IsotropicAntennaElement(‘FrequencyRange’,[1e9 2e9]);
fc = 1e9*(1:0.1:2);
resp = antenna(fc,[0;0]) Take the phased.IsotropoicAntennaElement as an example. ‘FrequencyRange’ is set between 1GHz and 2GHz. Under the default cofiguration, the response of the antenna in the specified ‘FrequencyRange’ is a constant of one. I wonder how can I change the respose with respect to different frequency.
antenna = phased.IsotropicAntennaElement(‘FrequencyRange’,[1e9 2e9]);
fc = 1e9*(1:0.1:2);
resp = antenna(fc,[0;0]) antenna response MATLAB Answers — New Questions
phased array antenna simulation
i want to simualte antenna radition pattern of radar rx antenna using 8 element linear antenna (as shown in figure 2).i could able to create the required radiation pattern by placing the antenna element on Y axis(ref figure 1) ,(Matlab code is attached for reference):
%%%%%%%%%START of matlab code for antenna element placing on Y axis %%%%%%%%%%%%%
clc;
clear all;
close all;
%fc = 77e9;
fc = 60e9;
%c = 3e8;
c = physconst(‘LightSpeed’);
lambda = c/fc;
Nt = 2;
Nr = 4;
ang = -90:90;
% % If both arrays have half-wavelength spacing, which are sometimes referred to as full arrays,
% % then the two-way pattern is close to the receive array pattern.
dr = lambda/2;
dt = Nr*(lambda/2);
Vy_pos1 = [0;-lambda;0];
Vy_pos2 = [0;-lambda+(lambda/2);0];
Vy_pos3 = [0;-lambda+ 2*(lambda/2);0];
Vy_pos4 = [0;-lambda+ 3*(lambda/2);0];
Vy_pos5 = [0;-lambda+ 4*(lambda/2);0];
Vy_pos6 = [0;-lambda+ 5*(lambda/2);0];
Vy_pos7 = [0;-lambda+ 6*(lambda/2);0];
Vy_pos8 = [0;-lambda+ 7*(lambda/2);0];
vxarray_4 = phased.ConformalArray(‘ElementPosition’,[Vy_pos1,Vy_pos2,Vy_pos3,Vy_pos4,Vy_pos5,Vy_pos6,Vy_pos7,Vy_pos8]);
figure
viewArray(vxarray_4 ,’ShowIndex’,’All’,’ShowNormal’,true,’Title’,’ conformal virtual array ,aantenna elements placed on Y axis’);
patv_4 = pattern(vxarray_4,fc,ang,0,’Type’,’powerdb’);
figure;
helperPlotMultipledBPattern(ang,[patv_4],[-30 0],…
{‘Two-way Pattern’,’Virtual Array Pattern’},…
‘Patterns of thin/full arrays and virtual array’,…
{‘-‘,’–‘},[1 2]);
%%%%%%%%%End of matlab code for antenna element placing on Y axis %%%%%%%%%%%%%
But my requirement is to achieve the same radaiation pattern by keeping antenna elements placement on X axis insetad of Y axis.please kindly check and confirm on how to realise same pattern by placing antenna element on the X axis(instead of y axis)?
i tried to place antenna elements on X axis and try to plot the resultant pattern but it doesnot meet the my required pattern.
figures and codes of the same are attached for reference.
%%%%%%%%%START of matlab code for antenna element placing on X axis %%%%%%%%%%%%%
clc;
clear all;
close all;
%fc = 77e9;
fc = 60e9;
%c = 3e8;
c = physconst(‘LightSpeed’);
lambda = c/fc;
Nt = 2;
Nr = 4;
ang = -90:90;
% % If both arrays have half-wavelength spacing, which are sometimes referred to as full arrays,
% % then the two-way pattern is close to the receive array pattern.
dr = lambda/2;
dt = Nr*(lambda/2);
Vx_pos1 = [-lambda;0;0];
Vx_pos2 = [-lambda+(lambda/2);0;0];
Vx_pos3 = [-lambda+ 2*(lambda/2);0;0];
Vx_pos4 = [-lambda+ 3*(lambda/2);0;0];
Vx_pos5 = [-lambda+ 4*(lambda/2);0;0];
Vx_pos6 = [-lambda+ 5*(lambda/2);0;0];
Vx_pos7 = [-lambda+ 6*(lambda/2);0;0];
Vx_pos8 = [-lambda+ 7*(lambda/2);0;0];
vxarray_3 = phased.ConformalArray(‘ElementPosition’,[Vx_pos1,Vx_pos2,Vx_pos3,Vx_pos4,Vx_pos5,Vx_pos6,Vx_pos7,Vx_pos8]);
figure
viewArray(vxarray_3 ,’ShowIndex’,’All’,’ShowNormal’,true,’Title’,’ conformal virtual array ,antenna elements placed on X axis’);
patv_3 = pattern(vxarray_3,fc,ang,0,’Type’,’powerdb’);
figure;
helperPlotMultipledBPattern(ang,[patv_3],[-30 0],…
{‘Two-way Pattern’,’Virtual Array Pattern’},…
‘Patterns of thin/full arrays and virtual array’,…
{‘-‘,’–‘},[1 2]);
%%%%%%%%%End of matlab code for antenna element placing on X axis %%%%%%%%%%%%%
please kindly check and confirm on how to realise rquired pattern by placing antenna element on the X axis(instead of y axis)?i want to simualte antenna radition pattern of radar rx antenna using 8 element linear antenna (as shown in figure 2).i could able to create the required radiation pattern by placing the antenna element on Y axis(ref figure 1) ,(Matlab code is attached for reference):
%%%%%%%%%START of matlab code for antenna element placing on Y axis %%%%%%%%%%%%%
clc;
clear all;
close all;
%fc = 77e9;
fc = 60e9;
%c = 3e8;
c = physconst(‘LightSpeed’);
lambda = c/fc;
Nt = 2;
Nr = 4;
ang = -90:90;
% % If both arrays have half-wavelength spacing, which are sometimes referred to as full arrays,
% % then the two-way pattern is close to the receive array pattern.
dr = lambda/2;
dt = Nr*(lambda/2);
Vy_pos1 = [0;-lambda;0];
Vy_pos2 = [0;-lambda+(lambda/2);0];
Vy_pos3 = [0;-lambda+ 2*(lambda/2);0];
Vy_pos4 = [0;-lambda+ 3*(lambda/2);0];
Vy_pos5 = [0;-lambda+ 4*(lambda/2);0];
Vy_pos6 = [0;-lambda+ 5*(lambda/2);0];
Vy_pos7 = [0;-lambda+ 6*(lambda/2);0];
Vy_pos8 = [0;-lambda+ 7*(lambda/2);0];
vxarray_4 = phased.ConformalArray(‘ElementPosition’,[Vy_pos1,Vy_pos2,Vy_pos3,Vy_pos4,Vy_pos5,Vy_pos6,Vy_pos7,Vy_pos8]);
figure
viewArray(vxarray_4 ,’ShowIndex’,’All’,’ShowNormal’,true,’Title’,’ conformal virtual array ,aantenna elements placed on Y axis’);
patv_4 = pattern(vxarray_4,fc,ang,0,’Type’,’powerdb’);
figure;
helperPlotMultipledBPattern(ang,[patv_4],[-30 0],…
{‘Two-way Pattern’,’Virtual Array Pattern’},…
‘Patterns of thin/full arrays and virtual array’,…
{‘-‘,’–‘},[1 2]);
%%%%%%%%%End of matlab code for antenna element placing on Y axis %%%%%%%%%%%%%
But my requirement is to achieve the same radaiation pattern by keeping antenna elements placement on X axis insetad of Y axis.please kindly check and confirm on how to realise same pattern by placing antenna element on the X axis(instead of y axis)?
i tried to place antenna elements on X axis and try to plot the resultant pattern but it doesnot meet the my required pattern.
figures and codes of the same are attached for reference.
%%%%%%%%%START of matlab code for antenna element placing on X axis %%%%%%%%%%%%%
clc;
clear all;
close all;
%fc = 77e9;
fc = 60e9;
%c = 3e8;
c = physconst(‘LightSpeed’);
lambda = c/fc;
Nt = 2;
Nr = 4;
ang = -90:90;
% % If both arrays have half-wavelength spacing, which are sometimes referred to as full arrays,
% % then the two-way pattern is close to the receive array pattern.
dr = lambda/2;
dt = Nr*(lambda/2);
Vx_pos1 = [-lambda;0;0];
Vx_pos2 = [-lambda+(lambda/2);0;0];
Vx_pos3 = [-lambda+ 2*(lambda/2);0;0];
Vx_pos4 = [-lambda+ 3*(lambda/2);0;0];
Vx_pos5 = [-lambda+ 4*(lambda/2);0;0];
Vx_pos6 = [-lambda+ 5*(lambda/2);0;0];
Vx_pos7 = [-lambda+ 6*(lambda/2);0;0];
Vx_pos8 = [-lambda+ 7*(lambda/2);0;0];
vxarray_3 = phased.ConformalArray(‘ElementPosition’,[Vx_pos1,Vx_pos2,Vx_pos3,Vx_pos4,Vx_pos5,Vx_pos6,Vx_pos7,Vx_pos8]);
figure
viewArray(vxarray_3 ,’ShowIndex’,’All’,’ShowNormal’,true,’Title’,’ conformal virtual array ,antenna elements placed on X axis’);
patv_3 = pattern(vxarray_3,fc,ang,0,’Type’,’powerdb’);
figure;
helperPlotMultipledBPattern(ang,[patv_3],[-30 0],…
{‘Two-way Pattern’,’Virtual Array Pattern’},…
‘Patterns of thin/full arrays and virtual array’,…
{‘-‘,’–‘},[1 2]);
%%%%%%%%%End of matlab code for antenna element placing on X axis %%%%%%%%%%%%%
please kindly check and confirm on how to realise rquired pattern by placing antenna element on the X axis(instead of y axis)? i want to simualte antenna radition pattern of radar rx antenna using 8 element linear antenna (as shown in figure 2).i could able to create the required radiation pattern by placing the antenna element on Y axis(ref figure 1) ,(Matlab code is attached for reference):
%%%%%%%%%START of matlab code for antenna element placing on Y axis %%%%%%%%%%%%%
clc;
clear all;
close all;
%fc = 77e9;
fc = 60e9;
%c = 3e8;
c = physconst(‘LightSpeed’);
lambda = c/fc;
Nt = 2;
Nr = 4;
ang = -90:90;
% % If both arrays have half-wavelength spacing, which are sometimes referred to as full arrays,
% % then the two-way pattern is close to the receive array pattern.
dr = lambda/2;
dt = Nr*(lambda/2);
Vy_pos1 = [0;-lambda;0];
Vy_pos2 = [0;-lambda+(lambda/2);0];
Vy_pos3 = [0;-lambda+ 2*(lambda/2);0];
Vy_pos4 = [0;-lambda+ 3*(lambda/2);0];
Vy_pos5 = [0;-lambda+ 4*(lambda/2);0];
Vy_pos6 = [0;-lambda+ 5*(lambda/2);0];
Vy_pos7 = [0;-lambda+ 6*(lambda/2);0];
Vy_pos8 = [0;-lambda+ 7*(lambda/2);0];
vxarray_4 = phased.ConformalArray(‘ElementPosition’,[Vy_pos1,Vy_pos2,Vy_pos3,Vy_pos4,Vy_pos5,Vy_pos6,Vy_pos7,Vy_pos8]);
figure
viewArray(vxarray_4 ,’ShowIndex’,’All’,’ShowNormal’,true,’Title’,’ conformal virtual array ,aantenna elements placed on Y axis’);
patv_4 = pattern(vxarray_4,fc,ang,0,’Type’,’powerdb’);
figure;
helperPlotMultipledBPattern(ang,[patv_4],[-30 0],…
{‘Two-way Pattern’,’Virtual Array Pattern’},…
‘Patterns of thin/full arrays and virtual array’,…
{‘-‘,’–‘},[1 2]);
%%%%%%%%%End of matlab code for antenna element placing on Y axis %%%%%%%%%%%%%
But my requirement is to achieve the same radaiation pattern by keeping antenna elements placement on X axis insetad of Y axis.please kindly check and confirm on how to realise same pattern by placing antenna element on the X axis(instead of y axis)?
i tried to place antenna elements on X axis and try to plot the resultant pattern but it doesnot meet the my required pattern.
figures and codes of the same are attached for reference.
%%%%%%%%%START of matlab code for antenna element placing on X axis %%%%%%%%%%%%%
clc;
clear all;
close all;
%fc = 77e9;
fc = 60e9;
%c = 3e8;
c = physconst(‘LightSpeed’);
lambda = c/fc;
Nt = 2;
Nr = 4;
ang = -90:90;
% % If both arrays have half-wavelength spacing, which are sometimes referred to as full arrays,
% % then the two-way pattern is close to the receive array pattern.
dr = lambda/2;
dt = Nr*(lambda/2);
Vx_pos1 = [-lambda;0;0];
Vx_pos2 = [-lambda+(lambda/2);0;0];
Vx_pos3 = [-lambda+ 2*(lambda/2);0;0];
Vx_pos4 = [-lambda+ 3*(lambda/2);0;0];
Vx_pos5 = [-lambda+ 4*(lambda/2);0;0];
Vx_pos6 = [-lambda+ 5*(lambda/2);0;0];
Vx_pos7 = [-lambda+ 6*(lambda/2);0;0];
Vx_pos8 = [-lambda+ 7*(lambda/2);0;0];
vxarray_3 = phased.ConformalArray(‘ElementPosition’,[Vx_pos1,Vx_pos2,Vx_pos3,Vx_pos4,Vx_pos5,Vx_pos6,Vx_pos7,Vx_pos8]);
figure
viewArray(vxarray_3 ,’ShowIndex’,’All’,’ShowNormal’,true,’Title’,’ conformal virtual array ,antenna elements placed on X axis’);
patv_3 = pattern(vxarray_3,fc,ang,0,’Type’,’powerdb’);
figure;
helperPlotMultipledBPattern(ang,[patv_3],[-30 0],…
{‘Two-way Pattern’,’Virtual Array Pattern’},…
‘Patterns of thin/full arrays and virtual array’,…
{‘-‘,’–‘},[1 2]);
%%%%%%%%%End of matlab code for antenna element placing on X axis %%%%%%%%%%%%%
please kindly check and confirm on how to realise rquired pattern by placing antenna element on the X axis(instead of y axis)? phased array antenna simulation MATLAB Answers — New Questions