Solving nonlinear equations using matrix

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!

[Torsten]

Specify all known values as numerical constants and try "fsolve" to solve your system.

I already did this for a bigger system that you posted - so you should know how to do it.

[Aquatris]

Take a look at this example here and you can adjust the code to your needs.

[Ismita]

Thanks @Aquatris

[Ismita]

Thanks @Torsten.

Answer 1

[Kartik Saxena]

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

[Ismita]

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