Solving equation in Matlab

Hi I obtained these results using Mathematica. If I want to get the same with Matlab, how should I do? At least I will get the red column Results :

The code :

Here is my code in Matlab :

clear all

L0= 941 * 10^9;

mi0=636 * 10^9;

k0= L0 + mi0 * (2/3);

kl = 900 * 10^9;

phi_m = 0.1;

Km = 10 ^-18;

w=2*pi;

eta=1000;

ac=10;

phi=0.12;

l=pi*0.001/phi;

t = (phi_m * (l^2)*eta)/(Km*kl);

L=(k-mi*(2/3));

gamma = 1 + (ac/pi) * ((k-k0)/(phi*mi)) * ((L+2*mi)/(L+mi));

For the rest :

that is : k=k0+ phi*k1 and µ=µ0 + phi*µ

[Kai Domhardt]

Have you had a look at the symbolic math toolbox?

[Pouyan Msgn]

I don't know how I can use it in this case. I am not professional and it is very complex

Answer 1

[Abraham Boayue]

See if you can make this code produce the results you want. The program by its self is right, but there might be some issues with the way I defined your equations and the values I chose for some constants.

% Fixed parameters
lamda_null  = 941e9;
  mue_null  = 636e9;
   phi_null = 0.1;
          w = 2*pi;
         k1 = 900e9;
         Km = 10 ^-18;
% Chosen parameters
eta  = [100 100 100 1000 1000 1000];
ac   =  10;
phi  =  [0.06 0.07 0.08 0.12 0.14 0.16];
phi_m = 0.1;
mue1  =   1;
N = length(eta);
Qp = zeros(1,N);
Qs = zeros(1,N);

% Calculated parameters
 k0     =  lamda_null + (2/3)*mue_null;
for i = 1:N
 kL     =  1e-6/phi(i);
 L      =  pi*0.001/phi(i);

   mue    =  mue_null + phi(i)*mue1;
   k      =  k0+ phi(i)*k1;
   tau    =  (phi_m *L^2*eta(i))/(Km*kL);

   % Define major equations 

   % 1. gamma
   lamda  =  (2/3)*(k-mue1);
   g1     =  (lamda+2*mue)/(lamda+mue);
   gamma  =  1+(ac/pi)*(k1/mue)*g1;

   % 2. k1/k0
   g2     =  (1+1i*w*tau)   / (1+1i*gamma*w*tau);
   g3     =  (1-(1-1/gamma) / (1+1i*w*tau));
   g4     =  ( 1 - ( 1i*kL^2 / gamma*w*tau)  * (1+1i*gamma*w*tau) );

   k1_k0  =  -4/3*(k0/mue)*g1*g2*(g3*g4).^(-1);

   % 3. mue1/mue0
   g5         = (lamda+2*mue)/(3*lamda+4*mue);
   g6         = (1+ 1i*4*ac*w*eta(i) / pi*mue);
   mue1_mue0  = -2/15*(8/3*g1*g2 +16*g5*(g5*g6).^(-1));

K  = k1_k0;
Me = mue1_mue0;
qp = imag(K+4/3*Me)/real(K+4/3*Me);
qs = imag(Me)/real(Me);
Qp(i) = qp;
Qs(i) = qs;
end
disp('Qp');
disp(Qp);
disp('Qs');
disp(Qs);

[Pouyan Msgn]

Thank you! The results I got were:

Qp 1.0e-14 *

   -0.6153   -0.6154   -0.6155   -0.0616   -0.0616   -0.0616

Qs 1.0e-14 *

   -0.6153   -0.6154   -0.6155   -0.0616   -0.0616   -0.0616

I have to work with this code. Thank you anyway