Solution of system of nonlinear equations

I have following system of equations to find q1, q2, q3 and q4, which I am unable to solve using solve using matlab 'solve' or 'fsolve' commands. All the values should be positive and should not be exceeding 1.0. Are there any suggestions? Or should I implement newton method to solve it?

eqc1_pre = 91.907*q1 - 1.1978e-14*q3 + 3.534e-14*q1*q3 + 9.1309e+6*(0.007649*q1 + 0.03937*q3)*(0.0038245*q1^2 + 0.03937*q1*q3 + 0.18186*q3^2 - 996.13*q3 - 0.03937*q4 + 1.6325e+7) + 1.819e-13*q3^2 == 0
eqc2_pre = 2914.4*q2 - 875.08*q3 + 3645.7*q2*q3 + 1.3054e+7*(0.007649*q2 + 0.03937*q3)*(0.0038245*q2^2 + 0.03937*q2*q3 + 0.18186*q3^2 - 996.23*q3 - 0.03937*q4 + 1.6325e+7) + 18765.0*q3^2 == 0
eqc3_pre = 3.638e-13*q1*q3 - 875.08*q2 - 9.4963e+8*q3 - 18765.0*q4 - 1.1978e-14*q1 + 37530.0*q2*q3 + 9.1309e+6*(0.03937*q1 + 0.36373*q3 - 996.13)*(0.0038245*q1^2 + 0.03937*q1*q3 + 0.18186*q3^2 - 996.13*q3 - 0.03937*q4 + 1.6325e+7) + 1.3054e+7*(0.03937*q2 + 0.36373*q3 - 996.23)*(0.0038245*q2^2 + 0.03937*q2*q3 + 0.18186*q3^2 - 996.23*q3 - 0.03937*q4 + 1.6325e+7) + 1.767e-14*q1^2 + 1822.8*q2^2 + 260050.0*q3^2 + 7.7808e+12 == 0
eqc4_pre = - 1374.8*q1^2 - 14153.0*q1*q3 - 1965.6*q2^2 - 20235.0*q2*q3 - 158850.0*q3^2 + 8.701e+8*q3 + 34387.0*q4 - 1.4259e+13 == 0
 

[Aquatris]

You can try numeric approaches like fmincon but I doubt you can solve this since you have values ranging from 1e-14 to 1e12, which many programing languages will have difficulty accuratly representing in a single equation.

[Pavl M.]

clc
clear all;

global st=0; %start of constraint intverval
global lst=1; % end of constraint interval

% 79% Completed:
function F = equations(q)
q1 = q(1);
q2 = q(2);
q3 = q(3);
q4 = q(4);
pi = 3.141592653589793;
F(1) = 91.907945087554152399491554310275*q1 - 6.9530743154589888866367962498927e-42*q3 + 0.000000000000035339972092069079645555751006121*q1*q3 + 0.0000000000001818989403545856475830078125*q3^2 + 9130872.0607999999462988460436463*((25*q1*pi^2)/32258 + (5*q3)/127)*((25*pi^2*q1^2)/64516 + (5*q1*q3)/127 + (2579522850115110167785284050996441184169*q3^2)/14183735516765013986772738621315219456000 - (269437024701681245919274115996940028669*q3)/270482702025760519408161296220160000 - (5*q4)/127 + 1010469987441913066378883825611918454307/61897001964269013744956211200000);
F(2) = 3061.5545338275630699055886034165*q2 - 0.00000000000000000000000052628096743060987461514329634383*q3 + 3645.7069886924359708335725693911*q2*q3 + 18764.877243218995339070902448522*q3^2 + 13054497.238043855951388792124012*((25*q2*pi^2)/32258 + (5*q3)/127)*((25*pi^2*q2^2)/64516 + (5*q2*q3)/127 + (2579522850115110167785284050996441184169*q3^2)/14183735516765013986772738621315219456000 - (269464072971883821971214932126562044669*q3)/270482702025760519408161296220160000 - (5*q4)/127 + 1010469987441913066378883825611918454307/61897001964269013744956211200000);
F(3) = 0.000000000000363797880709171295166015625*q1*q3 - 0.00000000000000000000000052628096743060987461514329634383*q2 - 949606892.52341601879283837341561*q3 - 18764.877243218995520969842803107*q4 - 1960*q3*(127/(10*pi) + 127/20) - 6.9530743154589888866367962498927e-42*q1 + 37529.754486437990678141804897043*q2*q3 + 9130872.0607999999462988460436463*((5*q1)/127 + (2579522850115110167785284050996441184169*q3)/7091867758382506993386369310657609728000 - 269437024701681245919274115996940028669/270482702025760519408161296220160000)*((25*pi^2*q1^2)/64516 + (5*q1*q3)/127 + (2579522850115110167785284050996441184169*q3^2)/14183735516765013986772738621315219456000 - (269437024701681245919274115996940028669*q3)/270482702025760519408161296220160000 - (5*q4)/127 + 1010469987441913066378883825611918454307/61897001964269013744956211200000) + 13054497.238043855951388792124012*((5*q2)/127 + (2579522850115110167785284050996441184169*q3)/7091867758382506993386369310657609728000 - 269464072971883821971214932126562044669/270482702025760519408161296220160000)*((25*pi^2*q2^2)/64516 + (5*q2*q3)/127 + (2579522850115110167785284050996441184169*q3^2)/14183735516765013986772738621315219456000 - (269464072971883821971214932126562044669*q3)/270482702025760519408161296220160000 - (5*q4)/127 + 1010469987441913066378883825611918454307/61897001964269013744956211200000) + 0.00000000000001766998604603453982277787550306*q1^2 + 1822.8534943462179854167862846955*q2^2 + 260045.56651039855910887591026306*q3^2 + 7780857152948.6118200978386072826;
F(4) = - 139.29999999999999918073863084391*pi^2*q1^2 - 199.15857468494441639131927665645*pi^2*q2^2 - 14152.879999999999916763044893742*q1*q3 - 20234.511187990352705358038508296*q2*q3 - 158847.98154891714318971301560184*q3^2 + 870095677.35924348389019942898346*q3 + 34387.391187990352622121083402037*q4 - 14258923878318.356547643662453739;
end

function q = steffensen_method(F, q0, tol, max_iter)
global st;
global lst;
q = q0;
for iter = 1:max_iter
    Fq = F(q);
    q1 = q - Fq;
    Fq1 = F(q1);
    q2 = q1 - Fq1;
    Fq2 = F(q2);
    q_new = q - (Fq.^2) ./ (Fq2 - 2*Fq1 + Fq);
    if norm(q_new - q) < tol && all (q_new > st & q_new <= lst)
        q = q_new;
        return;
    end
    q = q_new;
end
warning('Steffensen method did not converge');
end
function secant_method_optimized()
global st;
global lst;
% Initial guesses
q0 = [1; 1; 1; 1];  % First initial guess
q1 = [1.1; 1.1; 1.1; 1.1];  % Second initial guess
% Tolerance and maximum iterations
tol = 1e-5;
max_iter = 1500;
% Secant method iteration
for iter = 1:max_iter
    % Evaluate function at q0 and q1
    F0 = equations(q0);
    F1 = equations(q1);
    % Secant step
    q_next = q1 - (F1 .* (q1 - q0)) / (F1 - F0+eps);
    % Update guesses
    q0 = q1;
    q1 = q_next;
    % Check for convergence
    if norm(q1 - q0,1) < tol && all (q1 > st & q1 <= lst)
        fprintf('Converged to solution in %d iterations\n', iter);
        disp(q1);
        return;
    end
end
warning('Secant Method Optimized, Failed to converge in %d iterations\n', max_iter);
end

function q = fixed_point_iteration(q0, tol, max_iter)
global st;
global lst;
q = q0;
for iter = 1:max_iter
    q_new = equations(q);
    if norm(q_new - q,1) < tol && all (q_new > st & q_new <= lst)
        q = q_new;
        return;
    end
    q = q_new;
end
warning('Fixed-point iteration did not converge');
end
%To be continued to fix it:
function q = secant_method(q0, q1, tol, max_iter)
global st;
global lst;
for iter = 1:max_iter
    Fq0 = equations(q0);
    Fq1 = equations(q1);
    q2 = q1 - Fq1 * ((q1 - q0)./(Fq1 - Fq0+eps));
    q2(isinf(q2)) = 5000000;
    q2 = q2(:,1);
    %if sum(sqrt(abs(q2.^2 - q1.^2))) < tol % or norm(q2 - q1) < tol
    if norm(q2-q1,1) < tol && all (q_new > st & q_new <= lst)
        q = q2;
        fprintf('Converged in %d iterations.\n', iter);
        return;
    end
    q0 = q1;
    q1 = q2;
end
q = q2;
warning('Secant method did not converge');
end
function q = newton_m(initial_guess, tol, max_iter)

global st;
global lst;
q = initial_guess;
for iter = 1:max_iter
    F = equations(q);
    J = jacobian(q);
    delta_q = sum(-J \ F');
    q = q + delta_q';
    if norm(delta_q) < tol && all (q > st & q <= lst)
        fprintf('Converged in %d iterations.\n', iter);
        return;
    end
end
fprintf('Did not converge within %d iterations.\n', max_iter);
end

function q = broyden_method(initial_guess, tol, max_iter)

global st;
global lst;

q = initial_guess;
B = eye(length(q));  % Initial Jacobian approximation (identity matrix)

for iter = 1:max_iter
    F = equations(q);
    F = F(:);  % Ensure F is a column vector
    delta_q = -B \ F;  % Solve for delta_q
    q_new = q + delta_q;  % Update q
    % Check for convergence
    if norm(delta_q) < tol && all(q_new > st & q_new <= lst)
        fprintf('Converged in %d iterations.\n', iter);
        q = q_new;
        return;
    end
    F_new = equations(q_new);
    F_new = F_new(:);  % Ensure F_new is a column vector
    y = F_new - F;
    delta_q = delta_q(:); % Ensure delta_q is a column vector
    B = B + ((y - B * delta_q)'* delta_q) / (delta_q' * delta_q);  % Update Jacobian approximation
    q = q_new;  % Update q
end
fprintf('Did not converge within %d iterations.\n', max_iter);
end

function q = broyden_methodq(initial_guess, tol, max_iter)
global st;
global lst;
q = initial_guess;
B = eye(length(q));  % Initial Jacobian approximation (identity matrix)
for iter = 1:max_iter
    F = equations(q);
    delta_q = -B\F';
    q_new = q + sum(delta_q');
    if norm(delta_q) < tol && all (q_new > st & q_new <= lst)
        fprintf('Converged in %d iterations.\n', iter);
        return;
    end
    F_new = equations(q_new);
    y = F_new - F;
    B = B + ((y - B*delta_q).*delta_q)/(delta_q*delta_q');
    q = q_new;
end
fprintf('Did not converge within %d iterations.\n', max_iter);
end

function q = latin_pgd_method(initial_guess, tol, max_iter)
global st;
global lst;
% for LATIN-PGD method implementation
% This is a simplified version and may need adjustments for your specific problem
q = initial_guess
for iter = 1:max_iter
    % Solve the local problem
    q_local = solve_local_problem(q);
    % Solve the global problem
    q_global = solve_global_problem(q_local);
    % Update the solution
    q_new = update_solution(q, q_local, q_global);
    % Check for convergence
    if norm(q_new - q) < tol && all (q_new > st & q_new <= lst)
        fprintf('Converged in %d iterations.\n', iter);
        q = q_new;
        return;
        fprintf('LATIN-PGD method is not fully implemented in this example.\n');
    end
end
end
function q_local = solve_local_problem(q)
% Placeholder for solving the local problem
% Implement the local problem solution here
q_local = q;  % Example placeholder
end
function q_global = solve_global_problem(q_local)
% Placeholder for solving the global problem
% Implement the global problem solution here
q_global = q_local;  % Example placeholder
end
function q_new = update_solution(q, q_local, q_global)
% Placeholder for updating the solution
% Implement the solution update here
q_new = (q_local + q_global) / 2;  % Example placeholder
end
function u_corr = latin_pgd_complete(initial_guess)

global st;
global lst;
% Problem parameters
maxIter = 400;        % Maximum number of iterations
tol = 1e-6;           % Tolerance for convergence
% Initial guess
u = initial_guess;
% LATIN-PGD iteration
for iter = 1:maxIter
    % Prediction step
    u_pred = prediction_step(u);
    % Correction step
    u_corr = correction_step(u_pred);
    % Check for convergence
    if norm(u_corr - u,3) < tol && all (q_new > st & q_new <= lst)
        fprintf('Converged to solution in %d iterations\n', iter);
        #disp(u_corr);
        return;
    end
    % Update solution
    u = u_corr;
end
warning('Failed to converge in %d iterations\n', maxIter);
end
function u = initial_guess()
% Define an initial guess for the solution
u = rand(4, 1);  % Example with random initial guess, adapted to 4 variables
end
function u_pred = prediction_step(u)
% Prediction step of the LATIN method
% Example: using a simplified linear model for prediction
% Here we assume a linear relationship for illustration purposes
J = jacobian(u);
F = equations(u);
u_pred = u - J./(eps+F); %F*inv(J);  % Update with actual prediction logic
end
function u_corr = correction_step(u_pred)
% Correction step of the LATIN method
% Example: applying correction using nonlinear residuals
F = equations(u_pred);
u_corr = u_pred - 0.01 * F;  % Update with actual correction logic
end
function root = bisection_method(a, b, tol, max_iter,q)
global st;
global lst;
if equations(q+a) * equations(q+b)' >= 0
    error('The function must have different signs at a and b');
end
for iter = 1:max_iter
    c = (a + b) / 2;
    q = q+c;
    root = equations(q);
    if root == 0 || (b - a) / 2 < tol && all (q > st & q <= lst)
        fprintf('Converged in %d iterations.\n', iter);
        return;
    end
    if equations(q)*equations(q)' < 0
        b = c;
    else
        a = c;
    end
end
warning('Did not converge within the maximum number of iterations');
end
%Whether next Jacobians are correct (Fill in with your computed Jacobians):
function J = jacobian(q)
%{
 ##q1 = q(1);
 ##q2 = q(2);
 ##q3 = q(3);
 ##q4 = q(4);
%}
pi = 3.141592653589793;
J = zeros(4, 4);
n = length(q);
J = zeros(n, n);
h = 1e-6;  % Step size for finite differences
for i = 1:n
    for j = 1:n
        q1 = q;
        q2 = q;
        q1(j) = q1(j) + h;
        q2(j) = q2(j) - h;
        J(i, j) = (equations(q1)(i) - equations(q2)(i)) / (2 * h);
    end
end
% Partial derivatives for the first equation
%{
##J(1, 1) = 91.907945087554152399491554310275 + 0.000000000000035339972092069079645555751006121*q3 + 9130872.0607999999462988460436463*((25*pi^2)/32258)*((25*pi^2*q1^2)/64516 + (5*q1*q3)/127 + (2579522850115110167785284050996441184169*q3^2)/14183735516765013986772738621315219456000 - (269437024701681245919274115996940028669*q3)/270482702025760519408161296220160000 - (5*q4)/127 + 1010469987441913066378883825611918454307/61897001964269013744956211200000) + 9130872.0607999999462988460436463*((25*q1*pi^2)/32258)*((50*pi^2*q1)/64516 + (5*q3)/127); ##J(1, 2) = 0;
##J(1, 3) = -6.9530743154589888866367962498927e-42 + 0.000000000000035339972092069079645555751006121*q1 + 0.000000000000363797880709171295166015625*q3 + 9130872.0607999999462988460436463*((5)/127)*((25*pi^2*q1^2)/64516 + (5*q1*q3)/127 + (2579522850115110167785284050996441184169*q3^2)/14183735516765013986772738621315219456000 - (269437024701681245919274115996940028669*q3)/270482702025760519408161296220160000 - (5*q4)/127 + 1010469987441913066378883825611918454307/61897001964269013744956211200000) + 9130872.0607999999462988460436463*((25*q1*pi^2)/32258)*((5*q1)/127 + (5159045700230220335570568101992882368338*q3)/14183735516765013986772738621315219456000 - (269437024701681245919274115996940028669)/270482702025760519408161296220160000);
##J(1, 4) = -9130872.060799999946298846043646
##% Partial derivatives for the second equation
##J(2, 1) = 0;
##J(2, 2) = 3061.5545338275630699055886034165 + 3645.7069886924359708335725693911*q3 + 13054497.238043855951388792124012*((25*pi^2)/32258)*((25*pi^2*q2^2)/64516 + (5*q2*q3)/127 + (2579522850115110167785284050996441184169*q3^2)/14183735516765013986772738621315219456000 - (269464072971883821971214932126562044669*q3)/270482702025760519408161296220160000 - (5*q4)/127 + 1010469987441913066378883825611918454307/61897001964269013744956211200000) + 13054497.238043855951388792124012*((25*q2*pi^2)/32258)*((50*pi^2*q2)/64516 + (5*q3)/127);
##J(2, 3) = -0.00000000000000000000000052628096743060987461514329634383 + 3645.7069886924359708335725693911*q2 + 37529.754486437990678141804897043*q2 + 13054497.238043855951388792124012*((5)/127)*((25*pi^2*q2^2)/64516 + (5*q2*q3)/127 + (2579522850115110167785284050996441184169*q3^2)/14183735516765013986772738621315219456000 - (269464072971883821971214932126562044669*q3)/270482702025760519408161296220160000 - (5*q4)/127 + 1010469987441913066378883825611918454307/61897001964269013744956211200000) + 13054497.238043855951388792124012*((25*q2*pi^2)/32258)*((5*q2)/127 + (5159045700230220335570568101992882368338*q3)/14183735516765013986772738621315219456000 - (269464072971883821971214932126562044669)/270482702025760519408161296220160000);
##J(2, 4) = -13054497.238043855951388792124012*((25*q2*pi^2)/32258)*((5)/127);
##% Partial derivatives for the third equation
##%J(3, 1) = 0.000000000000363797880709171295166015625*q3 - 6.9530743154589888866367962498927e-42 + 9130872.0607999999462988460436463*((5)/127)*((25*pi^2*q1^2)/64516 + (5*q1*q3)/127 + (2579522850115110167785284050996441184169*q3^2)/14183735516765013986772738621315219456000 - (269437024701681245919274115996940028669*q3)/270482702025760519408161296220160000 - (5*q4)/127 + 1010469987441913066378883825611918454307/61897001964269013744956211200000) + 9130872.0607999999462988460436463*((25*q1*pi^2)/32258)*((5*q1)/127 + (5159045700230220335570568101992882368338*q3)/14183735516765013986772738621315219456000 - (269437024701681245919274115996940028669)/270482702025760519408161296220160000);
##%J(3, 2) = -0.00000000000000000000000052628096743060987461514329634383 + 37529.754486437990678141804897043*q3 + 1822.8534943462179854167862846955*q2;
##%J(3, 3) = 0.000000000000363797880709171295166015625*q1 - 949606892.52341601879283837341561 - 1960*(127/(10*pi) + 127/20) + 37529.754486437990678141804897043*q2 + 9130872.0607999999462988460436463*((5*q1)/127 + (5159045700230220335570568101992882368338*q3)/14183735516765013986772738621315219456000 - (269437024701681245919274115996940028669)/270482702025760519408161296220160000)*((25*pi^2*q1^2)/64516 + (5*q1*q3)/127 + (5159045700230220335570568101992882368338*q3)/14183735516765013986772738621315219456000 - (269437024701681245919274115996940028669)/270482702025760519408161296220160000 - (5*q4)/127 + 1010469987441913066378883825611918454307/61897001964269013744956211200000) + 13054497.238043855951388792124012*((5*q2)/127 + (5159045700230220335570568101992882368338*q3)/14183735516765013986772738621315219456000 - (269464072971883821971214932126562044669)/270482702025760519408161296220160000)*((25*pi^2*q2^2)/64516 + (5*q2*q3)/127 + (5159045700230220335570568101992882368338*q3)/14183735516765013986772738621315219456000
%Fill-in with your:
%J(3,4) =
%J(4,1) =
%J(4,2) =
%J(4,3) =
%J(4,4) =
%}
end


% Initial guess for the variables q1, q2, q3, q4
initial_guess = [0; 6; 1; -1];
% Solve the equations
%[q, fval, info] = fsolve(@equations, initial_guess);
tol = 1e-7;  % Tolerance for convergence
max_iter = 500;  % Maximum number of iterations
% Run the Newton's method
solution1 = newton_m(initial_guess, tol, max_iter);
disp('Solution Newton-Rhapson method:');
disp(solution1);
solution2 = latin_pgd_complete(initial_guess);
disp('Solution LATIN-PGD method:');
disp(solution2);
solution3 = broyden_method(initial_guess, tol, max_iter);
disp('Solution Broyden method:');
disp(solution3);
q0 = [0; 1; 1; 1];  % Initial guess
q1 = [2; 0; -2; 2];  % Second guess
tol = 1e-4;         % Tolerance
max_iter = 1000;     % Maximum number of iterations
solution4 = secant_method(q0, q1, tol, max_iter);
disp('Solution Secant method:')
disp(solution4);
q0 = [0; 1; -1; 2];  % Initial guess
tol = 1e-4;         % Tolerance
max_iter = 10000;     % Maximum number of iterations
solution5 = fixed_point_iteration(q0, tol, max_iter);
disp('Solution fixed point iteration methode:')
disp(solution5);
%solution_latin_pgd = latin_pgd_method(initial_guess, tol, max_iter);
%disp('Solution using LATIN-PGD method:');
%disp(solution_latin_pgd);
disp('Secant method optimized')
secant_method_optimized();
disp('Bisections methode:')
a = 1e-14;
b = 1e12;
root = bisection_method(a, b, tol, max_iter,q0);
disp(root)

%Constructed from needing help code by
%https://independent.academia.edu/PMazniker
%[phone]
%https://diag.net/u/u6r3ondjie0w0l8138bafm095b
%https://github.com/goodengineer
%https://orcid.org/0000-0001-8184-8166
%https://willwork781147312.wordpress.com/portfolio/cp/
%https://www.youtube.com/channel/UCC__7jMOAHak0MVkUFtmO-w
%https://nanohub.org/members/130066
%https://pangian.com/user/hiretoserve/
%https://substack.com/profile/191772642-paul-m

Please be attentive and careful with digital items, electronic artifacts and network logs and records.

If you correctly pay attention, review logs and records, you will see the truth.

I once completed full Newton Rhapson iterative procedure set for solving real world physics system of non-linear equations in TCE Matlab without any AI companionship and delivered it.

Who are they, others? What did they do? No, not that not. I do ok of course. Alternatives to Newton Rhapson:

There are several alternatives to the Newton-Raphson method for solving nonlinear equations.

I tested(checked) it in TCE.

A few specific functions to solution to your specific questions are fully completed by me. I need substantial support to complete it. It is ok and reasonable and even mission and normal life balancee supporting critical to hire me.

[John D'Errico]

@Pavl M.

Answers is not a job seeking forum. Please do not put requests for someone to hire you in your posts. I've removed them for you.

[Aquatris]

@Pavl M. Did you even solve the problem that was asked? The question askes for solutions within [0 1] space where as, if I understood your confusing comment, your solution does not put any constraints to the parameters.

Answer 1

[Alex Sha]

@Swami Between the range [0,1] you specified, most likely, there should be no real numbers solution，an approximate solution is:

q1: -31830.1002324947
q2: -31835.6758682725
q3: 6184.74933332136
q4: 336421377.881206

Answer 2

[John D'Errico]

This will be literally impossible to do in double precision, given the huge dynamic range of those coefficients. But we can look to see if it is even likely a solution does exist.

syms q [4,1]
eqc_pre(1) = 91.907*q1 - 1.1978e-14*q3 + 3.534e-14*q1*q3 + 9.1309e+6*(0.007649*q1 + 0.03937*q3)*(0.0038245*q1^2 + 0.03937*q1*q3 + 0.18186*q3^2 - 996.13*q3 - 0.03937*q4 + 1.6325e+7) + 1.819e-13*q3^2;
eqc_pre(2) = 2914.4*q2 - 875.08*q3 + 3645.7*q2*q3 + 1.3054e+7*(0.007649*q2 + 0.03937*q3)*(0.0038245*q2^2 + 0.03937*q2*q3 + 0.18186*q3^2 - 996.23*q3 - 0.03937*q4 + 1.6325e+7) + 18765.0*q3^2;
eqc_pre(3) = 3.638e-13*q1*q3 - 875.08*q2 - 9.4963e+8*q3 - 18765.0*q4 - 1.1978e-14*q1 + 37530.0*q2*q3 + 9.1309e+6*(0.03937*q1 + 0.36373*q3 - 996.13)*(0.0038245*q1^2 + 0.03937*q1*q3 + 0.18186*q3^2 - 996.13*q3 - 0.03937*q4 + 1.6325e+7) + 1.3054e+7*(0.03937*q2 + 0.36373*q3 - 996.23)*(0.0038245*q2^2 + 0.03937*q2*q3 + 0.18186*q3^2 - 996.23*q3 - 0.03937*q4 + 1.6325e+7) + 1.767e-14*q1^2 + 1822.8*q2^2 + 260050.0*q3^2 + 7.7808e+12;
eqc_pre(4) = - 1374.8*q1^2 - 14153.0*q1*q3 - 1965.6*q2^2 - 20235.0*q2*q3 - 158850.0*q3^2 + 8.701e+8*q3 + 34387.0*q4 - 1.4259e+13;

vpa(eqc_pre(:),4)

UGH. But we can look at whether it is likely a solution even exists in the hyper-box [0,1]^4.

double(subs(eqc_pre,q,rand(4,1)))

ans = 1×4

1.0e+17 * 0.0001 0.0001 -3.6066 -0.0001

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

double(subs(eqc_pre,q,rand(4,1)))

ans = 1×4

1.0e+17 * 0.0000 0.0000 -3.6070 -0.0001

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

double(subs(eqc_pre,q,rand(4,1)))

ans = 1×4

1.0e+17 * 0.0000 0.0000 -3.6077 -0.0001

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

However, I would point out that having substituted three sets of random numbers into those equations, almost always equation 3 tends to run about 3.6e17, and those computations were performed in high precision.

eq3 = matlabFunction(eqc_pre(3))

eq3 = function_handle with value:
    @(q1,q2,q3,q4)q1.*(-1.1978e-14)-q2.*8.7508e+2-q3.*9.4963e+8-q4.*1.8765e+4+q1.*q3.*3.638e-13+q2.*q3.*3.753e+4+(q1.*3.59483533e+5+q3.*3.321182257e+6-9.095563417e+9).*(q3.*(-9.9613e+2)-q4.*3.937e-2+q1.*q3.*3.937e-2+q1.^2.*3.8245e-3+q3.^2.*1.8186e-1+1.6325e+7)+q1.^2.*1.767e-14+q2.^2.*1.8228e+3+q3.^2.*2.6005e+5+(q2.*5.1393598e+5+q3.*4.74813142e+6-1.300478642e+10).*(q3.*(-9.9623e+2)-q4.*3.937e-2+q2.*q3.*3.937e-2+q2.^2.*3.8245e-3+q3.^2.*1.8186e-1+1.6325e+7)+7.7808e+12

eq3val = eq3(rand(1000000,1),rand(1000000,1),rand(1000000,1),rand(1000000,1));
min(eq3val)

ans = -3.6078e+17

max(eq3val)

ans = -3.6061e+17

This is enough to convince me that no solution exists in that region, since any of a set of 1e6 random points in the box NEVER see any deviation from a very narrow range, that is very far from zero.

Answer 3

[Pavl M.]

Here, I've solved it in more than 4 methods, for specific equations' system approximated solutions found with fixed interval constraints [a,b] added as per requirements plan:

A.1: Runs:

I've corrected, revamped functionality, amended the software program running codes presented in my comment above.

Working codes in Matlab compatible T.C.E. (fles SystemOfNonlinEqSolver1.m, newton_method.m, fronext.m, ...+):

Pictures are worth than thouthands words.I have the corresponding, according to specifications software program code maintaining running on file disk storage controller electronic unit processed machine. I developed them. They are for sale, if you need the software program codes running, contact me more + for more.