Solving nonlinear equations using matrix

1 次查看（过去 30 天）
Ismita
Ismita 2024-3-18
编辑： 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更早的评论

请先登录，再进行评论。

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

回答（1 个）

Kartik Saxena
Kartik Saxena 2024-3-21
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 个评论
Ismita
Ismita 2024-3-25
编辑：Ismita 2024-3-27
Ran in:
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

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Equation Solving 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by