Solving nonlinear equations using matrix

Question

Ismita 2024-3-18

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2095701-solving-nonlinear-equations-using-matrix

编辑： Ismita 2024-3-27

I have

p11 X1 + p12 X2 + p13 X3 + (-sin(psi) cos(X5) cos(X6) + cos(psi) sin(X5)) X4 &= C_1 (i) ;

p21 X1 + p22 X2 + p23 X3 + (-(1 + sin^2(psi)) cos(X5) cos(X6) + sin(psi) cos(psi) sin(X5)) X4 &= C_2(i) ;

p31 X1 + p32 X2 + p33 X3 + (-cos(X5) sin(X6)) X4 &= C_3 (i) \

p41 X1 + p42 X2 + p43 X3 + (- cos(psi) sin (psi) cos(X5) cos(X6) + (1 + cos^2(psi)) sin(X5)) X4 &= C_4(i) ;

p51 X1 + p52 X2 + p53 X3 + (- (1 + q) sin (psi) cos(X5)cos(X6) + (1 + q)cos(psi) sin(X5)) X4 &= C_5 (i) ;

p61 X1 + p62 X2 + p63 X3 - cos(X5) cos(X6) X4 &= C_6 (i);

Here i = 1: 6, C_1(i), C_2(i) ..... etc., coefficients p11, p12, ..... p63 are the constants and q = constant, psi = constant. Need to solve the equations for X1, X2, X3, X4, X5, X6. I am new to code to solve tis kind of equations. I request any suggestions and help. Thanks!

4 个评论
显示 2更早的评论隐藏 2更早的评论

Ismita 2024-3-18

Thanks @Aquatris

Ismita 2024-3-18

Thanks @Torsten.

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Kartik Saxena 2024-3-21

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2095701-solving-nonlinear-equations-using-matrix#answer_1428871

在 MATLAB Online 中打开

Hi,

Based on your problem, I'm adding a sample code snippet to solve this kind of problems, you can refer to it in order to design your solution according to your problem:

syms X1 X2 X3 X4 X5 X6 psi q
p = sym('p', [6, 3]); % Define the coefficient matrix p
C = sym('C', [6, 1]); % Define the constant vector C
eq1 = p(1,1)*X1 + p(1,2)*X2 + p(1,3)*X3 + (-sin(psi)*cos(X5)*cos(X6) + cos(psi)*sin(X5))*X4 == C(1);
eq2 = p(2,1)*X1 + p(2,2)*X2 + p(2,3)*X3 + (-(1 + sin(psi)^2)*cos(X5)*cos(X6) + sin(psi)*cos(psi)*sin(X5))*X4 == C(2);
eq3 = p(3,1)*X1 + p(3,2)*X2 + p(3,3)*X3 + (-cos(X5)*sin(X6))*X4 == C(3);
eq4 = p(4,1)*X1 + p(4,2)*X2 + p(4,3)*X3 + (-cos(psi)*sin(psi)*cos(X5)*cos(X6) + (1 + cos(psi)^2)*sin(X5))*X4 == C(4);
eq5 = p(5,1)*X1 + p(5,2)*X2 + p(5,3)*X3 + (-(1 + q)*sin(psi)*cos(X5)*cos(X6) + (1 + q)*cos(psi)*sin(X5))*X4 == C(5);
eq6 = p(6,1)*X1 + p(6,2)*X2 + p(6,3)*X3 - cos(X5)*cos(X6)*X4 == C(6);
% Solve the equations
sol = solve([eq1, eq2, eq3, eq4, eq5, eq6], [X1, X2, X3, X4, X5, X6]);
% Display the solutions
sol.X1
sol.X2
sol.X3
sol.X4
sol.X5
sol.X6

Also, refer to the following MathWorks documentations:

https://www.mathworks.com/help/symbolic/equation-solving.html?s_tid=CRUX_lftnav

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

Ismita 2024-3-25

编辑：Ismita 2024-3-27

在 MATLAB Online 中打开

NH_VOY2_Data.txt

thanks @Kartik Saxena. but it is not sloving.

clc;
clear; % all;
% Declare symbolic variables
syms X1 X2 X3 X4 X5 X6;
data = load('NH_VOY2_Data.txt');
% data = load('UxUyUz.txt');
mp = 1.67*10^(-24); %proton mass gm
mp_SI = 1.67*10^(-27); %proton mass kg
cc_m3 = 10^(-6); % m^3
k = 1.3807 * 10^(-16); % Boltzmann constant in CGS cm2 g s-2 K-1
% k = 1.380649*10^(-23); % boltzmann constant SI J/K
gamma = 5.00/3.00; % \gamma adiabatic index as 4/3 
time = data(:, 1);
density_CGS = data(:, 2); %number density
speed_km = data(:, 3); % 
temperature = data(:, 4); 
distance = data(:, 5); % upto 4 decimal
latitude = data(:, 6); % flow latitude angle theta 
longitude = data(:, 7); % flow longitude angle phi 
Ux = zeros(size(speed_km));
Uy = zeros(size(speed_km));
Uz = zeros(size(speed_km));
U = zeros(size(speed_km)); % magnitude of Ux, Uy, Uz
for i = 1:length(speed_km)    
    Ux(i) = speed_km(i) * cosd(latitude(i)) * cosd(longitude(i));
    Uy(i) = Ux(i) * tand(longitude(i));
    Uz(i) = speed_km(i) * sind(latitude(i));
    U(i) = sqrt(Ux(i)^2 + Uy(i)^2 + Uz(i)^2);       
end
%U_mag = sqrt(Ux.^2 + Uy.^2 + Uz.^2);
delta_umx = Ux - mean(Ux);
delta_umy = Uy - mean(Uy);
delta_umz = Uz - mean(Uz);
U_0_vec = mean(Ux) + mean(Uy) + mean(Uz); %equilibrium mean U_0 vector
U_0 = sqrt(mean(Ux)^2 + mean(Uy)^2 + mean(Uz)^2); %equilibrium mean U_0 magnitude
mass_density_CGS = density_CGS * mp; %\rho
rho_0 = mean(mass_density_CGS); %\rho_0 is mean mass density (at equilibrium)
delta_rho_m = mass_density_CGS - rho_0; %fluctuation in mass density
pressure_CGS = density_CGS .* k .* temperature; % P = nkT 
P_0 = mean(pressure_CGS);
delta_p_m = pressure_CGS - P_0;
a_0 = sqrt((gamma*P_0)/rho_0); % sound speed
epsilon_0 = (0.5 * U_0 * U_0) + a_0^2/(gamma - 1);
q = epsilon_0/U_0^2; 
M_0 = U_0/a_0;
psi = acos(mean(Uz)/U_0); % acos(x); % Calculate angle in radian by arccosine of x 
psi_degree = rad2deg(psi); % in degree
% Combine your data into a matrix where each column is a variable
data_matrix = [Ux, Uy, Uz]; % Each should be a column vector
% Number of observations (time points)
n = size(data_matrix, 1);
% Subtract the mean from each column (variable)
mean_subtracted_data = data_matrix - mean(data_matrix);
% Initialize the covariance matrix
cov_matrix = zeros(3, 3);
% Compute the covariance matrix manually
for i = 1:3
    for j = 1:3
        sum_product = sum(mean_subtracted_data(:, i) .* mean_subtracted_data(:, j));
        cov_matrix(i, j) = sum_product / (n - 1);
    end
end
% Display the manually computed covariance matrix
disp('Manually computed covariance matrix:');
Manually computed covariance matrix:
disp(cov_matrix)
    3.4330   -1.6459    0.1858
   -1.6459    6.5821   -2.2524
    0.1858   -2.2524   17.8136
% Continue from the previous code
% Perform eigenvalue decomposition on the covariance matrix
[eigenvectors, eigenvalues] = eig(cov_matrix);
% Extract the eigenvalues into a vector
eigenvalues_vector = diag(eigenvalues);
% Find the index of the smallest eigenvalue
[~, min_index] = min(eigenvalues_vector);
% The eigenvector corresponding to the smallest eigenvalue
min_variance_direction = eigenvectors(:, min_index);
% A = eigenvectors
% B = eigenvalues
% Display the eigenvector corresponding to minimum variance (this is orthogonal to the wave vector)
disp('Eigenvector corresponding to minimum variance:');
Eigenvector corresponding to minimum variance:
disp(min_variance_direction);
    0.9089
    0.4140
    0.0505
% Step 6: Compute wave vector components
kmx = min_variance_direction(1)
kmx = 0.9089
kmy = min_variance_direction(2)
kmy = 0.4140
kmz = min_variance_direction(3)
kmz = 0.0505
km = sqrt(kmx^2 + kmy^2 + kmz^2);
%a_0 in CGS km/s and wavevector km in km^-1
omega_m_prime = a_0 * km;
% omega_m_prime = a_0 * km/(10^5) ;
theta_a_plus = acos((1 / M_0) * (U_0* kmz) / omega_m_prime);
theta_a_minus = pi - theta_a_plus;
phi_a_plus = acos((U_0 * kmx) / (M_0 * sin(theta_a_plus) * omega_m_prime));
phi_a_minus = pi - phi_a_plus;
%equation of continuity
p_11 = 1
p_11 = 1
p_12 = 1 + (1/M_0)*sin(psi)*cos(phi_a_plus)*sin(theta_a_plus) + (1/M_0)*cos(psi)*cos(theta_a_plus)
p_12 = 1.4426e+03
p_13 = 1 - (1/M_0)*sin(psi)*cos(phi_a_minus)*sin(theta_a_minus) - (1/M_0)*cos(psi)*cos(theta_a_minus)
p_13 = 1.4426e+03
% Momentum $\hat{x}$ component,
p_21 = sin(psi)
p_21 = 0.9999
p_22 = sin(psi) + (1/M_0)*(1 + sin(psi)^2)*cos(phi_a_plus)*sin(theta_a_plus) + (1/M_0)*sin(psi)*cos(psi)*cos(theta_a_plus) + (1/M_0^2)*sin(psi)
p_22 = 2.5236e+06
p_23 = sin(psi) - (1/M_0)*(1 + sin(psi)^2)*cos(phi_a_minus)*sin(theta_a_minus) - (1/M_0)*sin(psi)*cos(psi)*cos(theta_a_minus) + (1/M_0^2)*sin(psi)
p_23 = 2.5236e+06
% Momentum $\hat{y}$ component,
p_31 = 0
p_31 = 0
p_32 = (1/M_0)*sin(phi_a_plus)*sin(theta_a_plus)
p_32 = 657.2918
p_33 = -(1/M_0)*sin(phi_a_minus)*sin(theta_a_minus)
p_33 = -657.2918
% Momentum $\hat{z}$ component,
p_41 = cos(psi)
p_41 = -0.0166
p_42 = cos(psi) + (1/M_0)*cos(psi)*sin(psi)*cos(phi_a_plus)*sin(theta_a_plus) + (1/M_0)*(1 + cos(psi)^2)*cos(theta_a_plus) + (1/M_0^2)*cos(psi)
p_42 = -4.1744e+04
p_43 = cos(psi) - (1/M_0)*cos(psi)*sin(psi)*cos(phi_a_minus)*sin(theta_a_minus) - (1/M_0)*(1 + cos(psi)^2)*cos(theta_a_minus) + (1/M_0^2)*cos(psi)
p_43 = -4.1744e+04
% Total energy equation
p_51 = 0.5
p_51 = 0.5000
p_52 = 0.5 + (1/M_0)*(1+ epsilon_0/U_0^2)*sin(psi)*cos(phi_a_plus)*sin(theta_a_plus) + (1/M_0)*(1+ epsilon_0/U_0^2)*cos(psi)*cos(theta_a_plus) + gamma/(M_0^2*(gamma - 1))
p_52 = 5.4577e+09
p_53 = 0.5 - (1/M_0)*(1+ epsilon_0/U_0^2)*sin(psi)*cos(phi_a_minus)*sin(theta_a_minus) - (1/M_0)*(1+ epsilon_0/U_0^2)*cos(psi)*cos(theta_a_minus) + gamma/(M_0^2*(gamma - 1))
p_53 = 5.4577e+09
% $\frac{\delta \hat{u}_{mx}}{U_0}$ equation gives, 
p_61 = 0;
p_62 = (1/M_0)*cos(phi_a_plus)*sin(theta_a_plus);
p_63 = -(1/M_0)*cos(phi_a_minus)*sin(theta_a_minus);
% p_64 = -cos(theta_v)*cos(phi_v);
% Given equation
%for i = 1:length(delta_rho_m)
C_1 = (delta_rho_m / rho_0) + sin(psi) * (delta_umx / U_0) + cos(psi) * (delta_umz / U_0);
C_2 = sin(psi) * (delta_rho_m / rho_0) + (1 + sin(psi)^2) * (delta_umx / U_0) + ...
      sin(psi) * cos(psi) * (delta_umz / U_0) + (sin(psi) / (M_0^2)) * (delta_p_m / (rho_0 * a_0^2));
C_3 = delta_umy / U_0;
% Given equation
C_4 = cos(psi) * (delta_rho_m / rho_0) + sin(psi) * cos(psi) * (delta_umx / U_0) + ...
      (1 + cos(psi)^2) * (delta_umz / U_0) + (cos(psi) / (M_0^2)) * (delta_p_m / (rho_0 * a_0^2));
C_5 = sin(psi) * (delta_umx / U_0) + cos(psi) * (delta_umz / U_0) + ...
      (a_0^2 / (U_0^2 * (gamma - 1))) * ((delta_p_m / P_0) - (delta_rho_m / rho_0)) + ...
      (epsilon_0 / (U_0^2)) * (delta_rho_m / rho_0 + sin(psi) * (delta_umx / U_0) + cos(psi) * (delta_umz / U_0));
C_6 = delta_umx / U_0;
% C_7 = delta_u_mz / U_0;
% C_8 = delta_rho_m / rho_0;
% C9 =  delta_p_m / (rho_0 * a_0^2);
%%%%%%%%%%%%%%%%%%%%%%
%syms X1 X2 X3 X4 X5 X6 psi q
%p = sym('p', [6, 3]); % Define the coefficient matrix p
%C = sym('C', [6, 1]); % Define the constant vector C
format long g;
p = [
    p_11, p_12, p_13;
    p_21, p_22, p_23;
    p_31, p_32, p_33;
    p_41, p_42, p_43;
    p_51, p_52, p_53;
    p_61, p_62, p_63]
p = 6x3
1.0e+00 *

                         1           1442.5783787177           1442.5783787177
         0.999862532251414          2523555.62438943          2523555.62438943
                         0           657.29178223393          -657.29178223393
       -0.0165806091501798         -41743.7318434778         -41743.7318434778
                       0.5          5457669447.96842          5457669447.96842
                         0          1443.10618410427          1443.10618410427
eq1 = p(1,1)*X1 + p(1,2)*X2 + p(1,3)*X3 + (-sin(psi)*cos(X5)*cos(X6) + cos(psi)*sin(X5))*X4 - C_1(1);
eq2 = p(2,1)*X1 + p(2,2)*X2 + p(2,3)*X3 + (-(1 + sin(psi)^2)*cos(X5)*cos(X6) + sin(psi)*cos(psi)*sin(X5))*X4 - C_2(1);
eq3 = p(3,1)*X1 + p(3,2)*X2 + p(3,3)*X3 + (-cos(X5)*sin(X6))*X4 - C_3(1);
eq4 = p(4,1)*X1 + p(4,2)*X2 + p(4,3)*X3 + (-cos(psi)*sin(psi)*cos(X5)*cos(X6) + (1 + cos(psi)^2)*sin(X5))*X4 - C_3(1);
eq5 = p(5,1)*X1 + p(5,2)*X2 + p(5,3)*X3 + (-(1 + q)*sin(psi)*cos(X5)*cos(X6) + (1 + q)*cos(psi)*sin(X5))*X4 - C_4(1);
eq6 = p(6,1)*X1 + p(6,2)*X2 + p(6,3)*X3 - cos(X5)*cos(X6)*X4 - C_5(1);
% Solve the equations
sol = solve([eq1, eq2, eq3, eq4, eq5, eq6], [X1, X2, X3, X4, X5, X6]);
% Display the solutions
sol.X1
 
ans =
 
Empty sym: 0-by-1
 
sol.X2
 
ans =
 
Empty sym: 0-by-1
 
sol.X3
 
ans =
 
Empty sym: 0-by-1
 
sol.X4
 
ans =
 
Empty sym: 0-by-1
 
sol.X5
 
ans =
 
Empty sym: 0-by-1
 
sol.X6
 
ans =
 
Empty sym: 0-by-1
 