System of nonlinear equations

Gencho · June 14th, 2017, 8:36 am

I am trying to repeat calculations from Hull(options futures and other derivatives) chapter "Using Equity Prices to Estimate Default Probabilities". I want to solve system of 2 non-linear equations with 2 unknowns:
\begin{cases}
E_0 = V_0 N(d_1) - L e^{-rt}N(d_2),\quad \quad V_0, \sigma_V -? \\
\sigma_E E_0 = N(d_1) \sigma_V V_0
\end{cases}
where $$d_2 = d_1 - \sigma_V \sqrt{t} = \frac{\log\frac{V_0}{K} + \frac{1}{2}\sigma^2t}{\sigma_V \sqrt{t}} - \sigma_V \sqrt{t}$$
But, using matlabs' fsolve function, even with feeding analytical Jacobian for different data inputs sometimes returns me:

No solutions found
or
Equation solved, fsolve stalled

This is even worse if I use more complex equations below. Could somebody give some hint how to solve this system?
or Could somebody give tips, ideas or reference how to choose good initial conditions for fsolve?
Finally, I eager to be able to solve more complex model than the written above.

Thanks!

This is actual problem I want to do:
$$
\begin{cases}
E_0 = V_0 N(x_1) - Le^{-rt} N(x_1 - \sigma \sqrt{t})
- V_0 (K/V_0)^{2\lambda} N(y_1) + Ke^{-rt}(K/V_0)^{2\lambda-2} N(y_1 - \sigma \sqrt{t})
\\
\sigma_E E_0 = \frac{\partial E_0}{ \partial V_0} \sigma_V V_0
\end{cases}
$$
The variables: $$x_1 = x_1(V_0,\sigma_V),y_1 = y_1(V_0,\sigma_V)$$
N(x) - cumulative distribution function of x, rest of variables are known.
I have added a plot of the sum of squares of the last system. From this figure, I assume that finding a solution is difficult because of quite a flat bottom.

https://ibb.co/hAK3v5[url=https://www.google.nl/url?sa=t&rct=j&q= ... ysymK0QGrA]https://ibb.co/hAK3v5[/url]

outrun · June 14th, 2017, 9:28 am

Just the other day I ran into Matlab's fsolve telling there is no solution, ..and it was right!
Have you checked that it can find a *known* solution?

You first equation looks like a Call with strike L and underlying V0? What the interpretation of the 2nd?

Gencho · June 14th, 2017, 9:36 am

Hi!
Thanks for super fast reply.
Exactly! First is call( in the second system it is a barrier call) to which we will apply ITO's lemma w.r.t. the value of companys' assets V. The 2nd equation flows from our assumption that value of the companys' equity follows geometric brownian motion:
\begin{equation}
dE_t = \mu_E E_t dt + \sigma_E E_t dW_t \label{eqn:1}
\end{equation}
This equation is useful, as instantaneous volatility of equity \sigma_E is observable from the market data. Estimation of volatility, is done by repeating the paragraph Estimating Volatility from Historical Data of Hull

E_t w.r.t. V (for notation convenience we have V_t=V) we obtain the following:
\begin{equation}
E_t = \frac{\partial E}{\partial V} dV + \frac{\partial^2 E}{\partial V^2} d \langle V \rangle \label{eqn:2}
\end{equation}
Calculating partial derivatives \frac{\partial E}{\partial V} = N(d_1), \frac{\partial^2 E}{\partial V^2} and after placing dV = \mu_V V dt + \sigma_V V dW_t and d \langle V \rangle = \sigma^2_V V^2 dt we pick the diffusion parts of above equations.We equalize them and derive:
\begin{equation}
\sigma_E E dW_t = N(d_1) \sigma_V V dW_t \label{eqn:3}
\end{equation}
$$\sigma_E E = N(d_1) \sigma_V V $$
This is described in Hull book, but I added some missing steps in the book of derivation. Ultimate goal: to do it for Barrier Option.

Oops, I did not answer your quesion. Yes, in case of first system of equations, everything works almost always. The problem appears with second system.

Cuchulainn · June 14th, 2017, 1:38 pm

Hull in a footnote in chapter 20.6 (6th ed.) mentions that this problem can be solved as a minimisation sum of squares. e.g. as in Excel Solve. It can also be solved I reckon by one of Matlab's DE solvers?

Is Matlab fsolve some kind of Newton thing? If the Jacobian becomes singular the problem breaks down.
2nd equation is Ito's lemma?

One possible outcome is that there _is_ a solution but that fsolve fails to converge for the particular set of parameters.Possible.

Cuchulainn · June 14th, 2017, 2:05 pm

Maybe you could post those numerical input values that you tested fsolve on?
To make things concrete.

kermittfrog · June 14th, 2017, 4:12 pm

Hi,

I remember doing stuff like this during uni, did it in Matlab as well, but do not remember the solver used...

Sensibile initial conditions:

Asset Level > Debt; some sensible leverage should suffice.
Asset volatility somewhat smaller than equity volatility.

Note that your model does not include any risk free rate or asset returns. Barrier options should be a simple extension.

Good luck.

Alan · June 14th, 2017, 7:10 pm

I gave it a shot in Mathematica, using Hull's example, and following his footnote suggestion to use a minimizer (in this case FindMinimum).

Worked fine.

The FindMinimum syntax is pretty self-evident, but the precise description is found here.
Note that E0 is the equity price while in Mathematica E is the base of the natural logs.

You can see my constraints and initial conditions in the code. (They were the first ones I tried, and mostly followed kermittfrog's suggestions).
The first element of the output is F^2 + G^2 at the found minimum:

I'm sure matlab must have something similar.

Clear[E0,sigE,r,T,L,Cumnormal,CN];
E0=3;
sigE=0.80;
r=0.05;
T=1;
L=10;
Cumnormal[xx_]:=(1+Erf[xx/Sqrt[2]])/2  
CN = Cumnormal; 

Clear[d1,d2,F,G];
d1[V0_,sigV_]:= (Log[V0/L] + (r + sigV^2/2)T)/(sigV Sqrt[T]);
d2[V0_,sigV_]:= d1[V0,sigV] - sigV Sqrt[T];
F[V0_,sigV_]:= E0 - V0 CN[d1[V0,sigV]] + L E^(-r T) CN[d2[V0,sigV]];
G[V0_,sigV_]:= CN[d1[V0,sigV]] sigV V0 - sigE E0;


FindMinimum[{F[V0,sigV]^2 + G[V0,sigV]^2,E0<V0,0<sigV<sigE},{{V0,5},{sigV,0.5 sigE}}]
{2.34285*10^-16,{V0->12.3954,sigV->0.212305}}

Cuchulainn · June 15th, 2017, 10:07 am

In C++11 with my hand-crafted DE solver

V = 12.3953871886
sigV = 0.212304713423

const double E0 = 3.0;
const double sigE = 0.80;
const double r = 0.05;
const double T = 1.0;
const double L = 10.0; // aka D == debt interest

// N(x) in C++11
auto CN = [](auto x) { return 0.5* (1.0 + std::erf(x / std::sqrt(2.0))); };

auto d1 = [&](auto V, auto sigV) 
 { return (std::log(V / L) + (r + sigV*sigV / 2)*T) / (sigV * std::sqrt(T));};
auto d2 = [&](auto V, auto sigV) { return d1(V, sigV) - sigV*std::sqrt(T); };

// F and G
auto F = [&](auto V, auto sigV)
{
 return E0 - V*CN(d1(V, sigV))+L*std::exp(-r*T)*CN(d2(V, sigV));
};

auto G = [&](auto V, auto sigV)
{
 return CN(d1(V, sigV))*sigV*V - sigE*E0;
};

double DefaultProbability(const ublas::vector<double>& trial) // 
{ // F^2 + G^2 

 double v = trial[0];
 double sig = trial[1];

 double f = F(v, sig);
 double g = G(v, sig);

 return f*f + g*g;
}

const int N_DIM = 2;
const int N_POP = 10000;
const int MAX_GENERATIONS = 10000;


int main()
{
 ublas::vector<double> solution;

 try 
 {
 ublas::vector<double> minBox(N_DIM);
 ublas::vector<double> maxBox(N_DIM);
 
 auto f = &DefaultProbability;
 DESolver<double> solver(N_DIM,N_POP, f);

 for (std::size_t i = 0; i < minBox.size(); ++i)
 {
 minBox[i] = 0.0;
 maxBox[i] = 20;
 }
 
 solver.Setup(minBox,maxBox,0.9,1.0);
 
 std::cout << "Calculating...\n\n";
 solver.Solve(MAX_GENERATIONS);
 solution = solver.Solution();

 std::cout << "\n\nBest Coefficients:\n";
 for (int i=0;i<solution.size();i++)
 {
 std::cout << i << ":" << std::setprecision(16) << solution[i] << std::endl;
 }

 
 }
 catch (std::exception& ex)
 {
 std::cout << ex.what();
 }

 

 return 0;

}

Cuchulainn · June 15th, 2017, 10:13 am

With an initial population = 20, and 100 generations we get

V = 12.39531192958

sigV = 0.2123048750404

One of the lessons is that iterative methods are too fragile for these kinds of problems (it is obvious if you have ever burned your fingers on NR). If there are mutiple solution they can converge to a local minimum at best. Differential Evolution always converges to a global minimum if such a minimum exists.

Things can break down in

1. Model
2. Numerical method
3. Programming errors

4. All three.

Gencho · June 15th, 2017, 11:28 am

Here is the slice of data I have problems with:
$$\sigma_E = 0,263$$
$$\text{st_debt} = 160911,467$$
$$\text{stock} = 33,846$$
$$\text{shares} = 1021,777$$
$$\text{lt_debt} = 153306,739$$
-----------------------------------------------------------------------------

Dear Cuchulainn,

- Exactly!
1) I was trying to find roots directly by matlab fsolve.
2) I have been following Hull suggestion( finding minimum of sum of squares of equations).

According to documentation: Matlab fsolve uses Levenberg-Marquardt algorithm which is basically a combination of gradient descent and Gauss-Newton methods. Comparing all algorithms Matlab has for fsolve, LMA provides most consistent results for my data.

I do not know if the Jacobian becomes singular( at least in command line I see no explicit information about that). But I can tell that I found Jacobian for the system, and I am feeding the solver with no luck( I see no difference except speed improvement).

Yes, 2nd equation is the Ito's lemma of E w.r.t. V.

I have been talking to supervisors( numerics is not their profile), they suggested to try 2 things:
1) to look not at minimization of f,g as:
$$
\min_x (y_1-f(x))^2 + (y_2-g(x))^2
$$
but at something like:
$$
\min_x (y_1-f(x))^4 + (y_2-g(x))^4
$$
or even higher powers.(I am checking this now)
2) use grid of 10x10 values of V, \sigma_V where possible solution might be situated and then compute the residual in each of these 100 points, and then start fsolve with the one that gave the smallest residual. (I have checked this, it did not work on data I post here.)
-----------------------------------------------------------------------------
Dear kermittfrog,

Yes, I use as initial condition for Asset level the following approach:
$$
V_0 = \sum \# \text{ outstanding stocks}_0\times \text{price of a stock}_0 + \text{short and long term debts}_0
$$
As initial value of the assets' volatility, \sigma_V, we use the half of the estimated above \sigma_E.

I think that risk free rate is included by r and yes, I simplify the model so that I have no asset returns.
-----------------------------------------------------------------------------

Dear Alan,

Thanks! My Matlab code also works well on Merton case. I am enquiring difficulties with the second system... To add, I have also tried feeding the solution of 1st system to the 2nd system.
-----------------------------------------------------------------------------

Is it meaningful to calculate condition number for this non-linear system?
Maybe there are some techniques to work with ill-posed problems?
I have been looking at the book by Heath ( Scientific Computing An Introductory Survey), where TOMS software is suggested as one of the strongest methods that implements homotopy/continuation methods which should always lead to global minimum.

So currently I will take a look at:
- Differential Evolution;
- find minimum of not sum of squares, but sum of higher powers.

P.S.I would like thanks everybody for posting ideas!

Cuchulainn · June 15th, 2017, 11:46 am

LevMar, that's another kettle of fish. It's awful. meshuggah.

You're welcome. Both Alan and myself have solved this (benign) problem for you, yes? QED.

find minimum of not sum of squares, but sum of higher powers.
Can't see why. Statisticians might do this (don't know why) but not really in mathematics. IMO it would be wasting time.

Is it meaningful to calculate condition number for this non-linear system?
No. BTW nonlinear systems don't have condition numbers AFAIK. Maybe you are thinking of matrices or sumthing.

"Maybe there are some techniques to work with ill-posed problems?"
There are but ...
It is not an ill-posed problem/ fsolve is the problem. IMO.

Cuchulainn · June 15th, 2017, 12:12 pm

To add, I have also tried feeding the solution of 1st system to the 2nd system.
What do you mean by this? Sounds scary..

Gencho · June 15th, 2017, 12:12 pm

Dear Cuchulainn,

Could you explain what do you mean by LevMar comment?

Yes, both of you have solved Merton case( 1st system of equations), which I also can solve. Problems appear with the second one.

Hehe, you are right. Statistician advised me on raising to higher powers.

Haha, I mean that the first system of equations is solvable. I was plugging output values for V, \sigma_V from the 1st system to the initial conditions of the 2nd system. Is it clear?

Here is my code, if someone wants to see failure of fsolve (code should be adjusted, as you have all values in the code of cell type, so switch cell to double):

global K H e_00 sigma_ee r s

sigma_e = sigma_e_6M{1}; 
V_0_prime_vector = st_debt{1} + lt_debt{1} + shares{1}.*stock{1};
sigma_v_prime_vector = sigma_e/2;
D = st_debt{1}+ 0.5*lt_debt{1};
s = 1;
r=0.05;
k=0;
gamma = 0;
e_0 = shares{1}.*stock{1};
list_DP=NaN(length(st_debt{1}),1);
list_SigmaV = NaN(length(st_debt{1}),1);
list_V = NaN(length(st_debt{1}),1);

options = optimoptions('fsolve');
options.Algorithm ='Levenberg-Marquardt';
options.TolX = 1e-6;
options.TolFun = 1e-6;
options.MaxIter = 10000;
options.MaxFunEvals = 10000;
options.Jacobian = 'on';
start_time = 1;
end_time = 1;
for i=start_time:end_time
        Dd = D(i);
        K = Dd;
        H = Dd;
        e_00 = e_0(i);
        sigma_ee = sigma_e(i);
        x0 =[V_0_prime_vector(i),sigma_e(i)] ;
        x = fsolve(@root2dBCrut,(x0),options);
        %--------------------------------------------------------
        nu = r - k - gamma - x(2) * x(2)/2;
        x_1 = ( log(x(1)/K) + nu * s )/(x(2) * sqrt(s));
        a_tilde = ( nu - gamma ) /(x(2)*x(2));
        x_2 = (H/x(1))^(2*a_tilde);
        x_3 = ( 2*log(H) - log(K*x(1)) + nu * s )/(x(2) * sqrt(s));
        sp = normcdf(x_1) - x_2 * normcdf(x_3);
        %---------------------------------------------------------
        list_DP(i) = 1-sp;
 end

where system of equations is:

function [F, Jac]  = root2dBCrut(x)
global K H e_00 sigma_ee r
lambda = 1 + r/(x(2)*x(2));
ers = exp( - r * s );
z_1 = ( log(x(1)/K) + lambda * x(2) * x(2) * s )/(x(2)*sqrt(s));
z_12 = z_1 - x(2) * sqrt(s);
z_2 = ( log(H*H/(x(1)*K)) + lambda * x(2) * x(2) * s )/(x(2)*sqrt(s));
z_22 = z_2 - x(2) * sqrt(s);

cdf_z_1 = normcdf(z_1);
cdf_z_12 = normcdf(z_12);
cdf_z_2 = normcdf(z_2);
cdf_z_22 = normcdf(z_22);

z_1_prime_V = 1/(x(1) * x(2) * sqrt(s));
z_12_prime_V = z_1_prime_V;
z_2_prime_V = z_1_prime_V;
z_22_prime_V = z_1_prime_V;

cdf_z_1_prime_V = normpdf(z_1) * z_1_prime_V;
cdf_z_12_prime_V = normpdf(z_12) * z_12_prime_V;
cdf_z_2_prime_V = normpdf(z_2) * z_2_prime_V;
cdf_z_22_prime_V = normpdf(z_22) * z_22_prime_V;

E = x(1) * cdf_z_1 - K * ers * cdf_z_12 - x(1) * (H/x(1))^(2*lambda) * cdf_z_2 + K * ers * (H/x(1))^(2*lambda-2) * cdf_z_22 ; 

a = cdf_z_1 + x(1) * cdf_z_1_prime_V;
b = - K * ers * cdf_z_12_prime_V;
c = - (H/x(1))^(2*lambda) * ( cdf_z_2 * (1-2*lambda) + x(1) *  cdf_z_2_prime_V );
d = K * ers * (H/x(1))^(2*lambda-2) * ( cdf_z_2 * (2-2*lambda)/x(1) + cdf_z_22_prime_V );

dEdV = a + b + c + d;

F(1) = e_00   - E;
F(2) = sigma_ee * e_00 -  x(1) * x(2) * dEdV;

Jac = SuperJacobianV3(x(1),x(2));
end

and Jacobian is:

function out = SuperJacobianV3(V,sV)
global Dd
K = Dd;
H =Dd;
r=0.05;
s=1;
lambda = ( r/2 + sV * sV/2 )/( sV * sV);
%------------------------------------------------------
x_1 = (  log(V/H) + (lambda * sV * sV * s )  ) / ( sV * sqrt(s) );
y_1 = ( -log(V/H) + (lambda * sV * sV * s )  ) / ( sV * sqrt(s) );
x_12 = x_1-sV* sqrt(s);
y_12 = y_1-sV* sqrt(s);
%------------------------------------------------------
x_1_prime_V = 1/(V * sV * sqrt(s));
y_1_prime_V = - 1/(V * sV * sqrt(s));
x_12_prime_V = x_1_prime_V;
y_12_prime_V = y_1_prime_V;
%------------------------------------------------------
x_1_prime_sV = - log(V/H)/ ( sqrt(s) * sV * sV ) - r * sqrt(s) / ( sV * sV ) + sqrt(s)/2;
y_1_prime_sV = - log(H/V) /( sqrt(s) * sV * sV ) - r * sqrt(s) / ( sV * sV ) + sqrt(s)/2;
x_12_prime_sV = - log(V/H)/ ( sqrt(s) * sV * sV ) - r * sqrt(s) / ( sV * sV ) - sqrt(s)/2;
y_12_prime_sV = - log(H/V) /( sqrt(s) * sV * sV ) - r * sqrt(s) / ( sV * sV ) - sqrt(s)/2;
%------------------------------------------------------
ers = exp( - r * s );
%------------------------------------------------------
cdf_x_1 = normcdf(x_1);
cdf_y_1 = normcdf(y_1);
cdf_x_12 = normcdf(x_12);
cdf_y_12 = normcdf(y_12);
%------------------------------------------------------
cdf_x_1_prime_V =    normpdf(x_1) * x_1_prime_V;
cdf_y_1_prime_V =    normpdf(y_1) * y_1_prime_V;
cdf_x_12_prime_V =   normpdf(x_12) * x_12_prime_V;
cdf_y_12_prime_V =   normpdf(y_12) * y_12_prime_V;
cdf_x_12_prime_V_prime_V = normpdf(x_12) * -x_12 * x_12_prime_V * x_12_prime_V + normpdf(x_12) * - x_1_prime_V / V;
%------------------------------------------------------
cdf_x_1_prime_sV = normpdf(x_1) * x_1_prime_sV;
cdf_y_1_prime_sV = normpdf(y_1) * y_1_prime_sV;
cdf_x_12_prime_sV = normpdf(x_12) * x_12_prime_sV;
cdf_y_12_prime_sV = normpdf(y_12) * y_12_prime_sV;
%--------------------------------------------------------------------------
dEdVp1 = cdf_x_1;
dEdVp2 = V * cdf_x_1_prime_V;
dEdVp3 = - K * ers * cdf_x_12_prime_V;
dEdVp4 = - (H/V)^(2*lambda) * V * ( 0 + (1-2*lambda)/V) * cdf_y_1;% assume mistake is here. if it is so, then it might be below symmetrically and in dvdsv and dvdv terms as well
dEdVp5 = - (H/V)^(2*lambda) * V * cdf_y_1_prime_V; 
dEdVp6 = K * ers * (H/V)^(2*lambda-2) * ( 0 + (2-2*lambda)/V) * cdf_y_12;
dEdVp7 = K * ers * (H/V)^(2*lambda-2) * cdf_y_12_prime_V;
dEdV = dEdVp1 + dEdVp2 + dEdVp3 + dEdVp4 + dEdVp5 + dEdVp6 + dEdVp7;
e11 = - dEdV;
%--------------------------------------------------------------------------
dEdsVp1 = V * cdf_x_1_prime_sV;
dEdsVp2 = - K * ers * cdf_x_12_prime_sV;
dEdsVp3 = - (H/V)^(2*lambda) * V * ( log(H/V) * 2 * - 2 * r/(sV^3) ) * cdf_y_1;
dEdsVp4 = - (H/V)^(2*lambda) * V * cdf_y_1_prime_sV; 
dEdsVp5 = K * ers * (H/V)^(2*lambda-2) * ( log(H/V) * 2 * - 2 * r/(sV^3)) * cdf_y_12;
dEdsVp6 = K * ers * (H/V)^(2*lambda-2) * cdf_y_12_prime_sV;
dEdsV = dEdsVp1 + dEdsVp2 + dEdsVp3 + dEdsVp4 + dEdsVp5 + dEdsVp6;
e12 = - dEdsV;
%--------------------------------------------------------------------------
ddEEdVVp1 = cdf_x_1_prime_V;
ddEEdVVp2 = cdf_x_1_prime_V - V * normpdf(x_1) * x_1_prime_V * ( x_1 * x_1_prime_V + 1/V );
ddEEdVVp3 = - K * ers * cdf_x_12_prime_V_prime_V;
ddEEdVVp4 = - (1-2*lambda) * (H/V)^(2*lambda) * ( cdf_y_1_prime_V - 2 * lambda * cdf_y_1/V);
ddEEdVVp5 = - (H/V)^(2*lambda) * V * ( ( 1-2*lambda)/V * cdf_y_1_prime_V  - normpdf(y_1) * y_1_prime_V * ( y_1 * y_1_prime_V + 1/V )  );% till this point above is checked
ddEEdVVp6 = K * ers * (H/V)^(2*lambda-2) * (2-2*lambda) * ( (1-2*lambda)*cdf_y_12 /(V*V) + cdf_y_12_prime_V/V );% this is  checked
ddEEdVVp7 = K * ers * (H/V)^(2*lambda-2) * ( (2-2*lambda)/V * cdf_y_12_prime_V - normpdf(y_12) * y_12_prime_V * ( y_12 * y_12_prime_V -1/V ) );% this is  checked
ddEdVV = ddEEdVVp1 + ddEEdVVp2 + ddEEdVVp3 + ddEEdVVp4 + ddEEdVVp5 + ddEEdVVp6 + ddEEdVVp7;
e21A1 = - sV * dEdV;
e21A2 = - sV * V * ddEdVV;
e21 = e21A1 + e21A2;
%--------------------------------------------------------------------------
ddEEdVsVp1 = cdf_x_1_prime_sV; 
ddEEdVsVp2 = - V * normpdf(x_1) * x_1_prime_V * (  x_1 * x_1_prime_sV + 1/sV ); % ok
ddEEdVsVp3 = K * ers * normpdf(x_12) * x_12_prime_V * (  x_12 * x_12_prime_sV + 1/sV ); % ok
ddEEdVsVp4 = (H/V)^(2*lambda) * ( - cdf_y_1 * 4 * r * ( log(H/V) * (2*lambda-1) + 1 )/(sV*sV*sV) + (2*lambda-1) * normpdf(y_1) * y_1_prime_sV );
ddEEdVsVp5 = V * (H/V)^(2*lambda) *  normpdf(y_1) * y_1_prime_V * ( y_1 * y_1_prime_sV + 1/sV + log(H/V) * 4 * r/(sV*sV*sV) );
ddEEdVsVp6 = K * ers * (H/V)^(2*lambda-2)/V * ( cdf_y_12 * 4 * r/(sV*sV*sV) * ( 1 + (2*lambda-2) * log(H/V)) + (2-2*lambda) * normpdf(y_12) * y_12_prime_sV );
ddEEdVsVp7 = - K * ers * (H/V)^(2*lambda-2) * normpdf(y_12) * y_12_prime_V * (  y_12 * y_12_prime_sV + 1/sV + 4 * r * log(H/V)/(sV*sV*sV) );
ddEdVsigmaV = ddEEdVsVp1 + ddEEdVsVp2 + ddEEdVsVp3 + ddEEdVsVp4 + ddEEdVsVp5 + ddEEdVsVp6 + ddEEdVsVp7;
e22A1 = - V * dEdV;
e22A2 = - sV * V * ddEdVsigmaV;
e22 = e22A1 + e22A2;
%--------------------------------------------------------------------------
out = [[e11 e12]; [e21 e22]];

Cuchulainn · June 15th, 2017, 12:21 pm

Could you post the Second system of equations? Did I miss something? Or have the goalposts been moved?

Cuchulainn · June 15th, 2017, 12:21 pm

Hehe, you are right. Statistician advised me on raising to higher powers.

LOL

Show them this

https://en.wikipedia.org/wiki/Runge%27s_phenomenon

System of nonlinear equations

System of nonlinear equations

Re: System of nonlinear equations

Re: System of nonlinear equations

Re: System of nonlinear equations

Re: System of nonlinear equations

Re: System of nonlinear equations

Re: System of nonlinear equations

Re: System of nonlinear equations

Re: System of nonlinear equations

Re: System of nonlinear equations

Re: System of nonlinear equations

Re: System of nonlinear equations

Re: System of nonlinear equations

Re: System of nonlinear equations

Re: System of nonlinear equations