How to Implement 3rd order equation in matlab script?

I want to Implement this equation in matlab but I am not getting how to do it. Can anyone help me with this? psi1,psi2 and K1 are dependent on theta.

[Sam Chak]

If both θ and Ψ are independent variables, then i is a function of two variables that produces a surface.

For example:

L                  = 160*membrane(1,100);
f                  = figure;
ax                 = axes;
s                  = surface(L); grid on
s.EdgeColor        = 'none';
view(3)
ax.XLim            = [1 201];
ax.YLim            = [1 201];
ax.ZLim            = [-53.4 160];
ax.CameraPosition  = [-145.5 -229.7 283.6];
ax.CameraTarget    = [77.4 60.2 63.9];
ax.CameraUpVector  = [0 0 1];
ax.CameraViewAngle = 36.7;
l1                 = light;
l1.Position        = [160 400 80];
l1.Style           = 'local';
l1.Color           = [0 0.8 0.8];
l2                 = light;
l2.Position        = [.5 -1 .4];
l2.Color           = [0.8 0.8 0];
s.FaceColor        = [0.9 0.2 0.2];
xlabel('\theta', 'fontsize', 16), ylabel('\Psi', 'fontsize', 16), zlabel('i', 'fontsize', 16)

Answer 1

[Jon]

I would do it something like this:

function ival = ifun(theta,psi)
% Define constants
K2 = 25; % just for example, you put in actual value
K3 = 32.245; % just for example, you put in actual value
% Determine values of constants m and n based on values of psi
if psi > psi1(theta)
    m = 1;
else
    m = 0;
end
if psi > psi2(theta);
    n = 1;
else
    n = 0;
end
% Evaluate overall function
ival = K1(theta)*psi + m*K2*(psi - psi1(theta))^2 + n*K3*(psi - psi2(theta));

% Define helper functions
    function k1val = K1(theta)
        k1val = theta^2 + 3; % just for example you put in real function here
    end
    function psi1Val = psi1(theta)
        psi1Val = 3*theta + theta^3; % just for example you put in real function here
    end
    function psi2Val = psi2(theta)
        psi2Val = exp(theta) + sin(theta);% just for example you put in real function here
    end
end

Then to actually evaluate ifun, either in a script or another function you would call it, for example

theta = pi/3;
psi = 0.9;
ival = ifun(theta,psi)

[Jon]

Note, I wouldn't use "i" as a variable name, as usually i is considered just a local array index, not a meaningful variable name, also i can be confused as indicating an imaginary number

[Prajwal Magadi]

Hello Jon,

Sorry for incomplete info, Actually psi1,psi2 and k1 are arrays dependent on theta. and K2=11,K3=185. after ivaluating I want to plot i vs psi.

Thank You

[Jon]

You could do it like this:

% Plot i as function of psi for fixed value of theta
theta = 20; % put your value of interest here
psi = linspace(0,1); % put your vector of psi values here

% loop to evaluate i at each value of psi
ivals = zeros(numel(psi),1); % preallocate
for k = 1:numel(psi)
    ivals(k) = ifun(theta,psi(k));
end

% plot results
plot(psi,ivals)
xlabel('psi')
ylabel('i')
title(['i(theta,psi) for theta = ',num2str(theta)])

% Define function to evaluate i, could put in seperate m file, if you want to reuse
% it for other applications
function ival = ifun(theta,psi)
% Assign constants
K2 = 11; % just for example, you put in actual value
K3 = 185; % just for example, you put in actual value

% Assign theta dependent parameters (or perhaps read in from a data file if it can
% change)
thetaVals = 0:3:30;
K1Vals = [67,62.5,53.5,38,23.5,17,14,12,10,8.75,8];
psi1Vals = [0.25,0.25,0.25,0.175,0.2,0.225,0.335,0.46,0.47,0.485,0.485];
psi2Vals = [0.25,0.25,0.25,0.25,0.275,0.35,0.43,0.495,0.545,0.56,0.56];

% Look up values of theta dependent parameters
K1 = interp1(thetaVals,K1Vals,theta);
psi1 = interp1(thetaVals,psi1Vals,theta);
psi2 = interp1(thetaVals,psi2Vals,theta);

% Determine values of constants m and n based on values of psi
if psi > psi1
    m = 1;
else
    m = 0;
end
if psi > psi2
    n = 1;
else
    n = 0;
end
% Evaluate overall function
ival = K1*psi + m*K2*(psi - psi1)^2 + n*K3*(psi - psi2);
end

I see that you have theta in degrees. Usually physical equations are defined in terms of radians. You should check whether theta and psi should be in radians or degrees, and make sure that you are consistent.

[Prajwal Magadi]

Thanks for your reply Jon. Here you have varied 'psi' values to get 'i' values. But I wants them in other way , vary 'i' values and get 'psi' values. And 'i'=1:18.

[Jon]

Here is how you could solve for psi given i and plot results. If you don't actually have to solve for exact values of psi and just want the data plotted the other way you could just modigy the code I provided previously, swapping the arguments to plot and the xlabel and ylabel, to plot psi vs i instead of i vs psi.

% Plot psi as function of i for fixed value of theta
theta = 20; % put your value of interest here
ivals = 1:18; % range of i to be plotted



% loop to evaluate psi at each value of i
psi = zeros(numel(ivals),1); % preallocate
for k = 1:numel(ivals)
    ival = ivals(k);
    % Define function whose roots will give you values of psi
    % need to do it inside of loop so that ival will be in scope
    zfun = @(psi) ival - ifun(theta,psi);
    psi(k) = fzero(zfun,[0,1]);
end

% plot results
figure
plot(ivals,psi)
ylabel('psi')
xlabel('i')
title(['psi(theta,i) for theta = ',num2str(theta)])

% Define function to evaluate i, could put in separate m file, if you want to reuse
% it for other applications
function ival = ifun(theta,psi)
% Assign constants
K2 = 11; % just for example, you put in actual value
K3 = 185; % just for example, you put in actual value

% Assign theta dependent parameters (or perhaps read in from a data file if it can
% change)
thetaVals = 0:3:30;
K1Vals = [67,62.5,53.5,38,23.5,17,14,12,10,8.75,8];
psi1Vals = [0.25,0.25,0.25,0.175,0.2,0.225,0.335,0.46,0.47,0.485,0.485];
psi2Vals = [0.25,0.25,0.25,0.25,0.275,0.35,0.43,0.495,0.545,0.56,0.56];

% Look up values of theta dependent parameters
K1 = interp1(thetaVals,K1Vals,theta);
psi1 = interp1(thetaVals,psi1Vals,theta);
psi2 = interp1(thetaVals,psi2Vals,theta);

% Determine values of constants m and n based on values of psi
if psi > psi1
    m = 1;
else
    m = 0;
end
if psi > psi2
    n = 1;
else
    n = 0;
end
% Evaluate overall function
ival = K1*psi + m*K2*(psi - psi1)^2 + n*K3*(psi - psi2);
end

[Jon]

Did this answer your question? If so please accept the answer so others will know a solution is available. Otherwise what remaining issues do you have?

Answer 2

[Sam Chak]

Hi @Prajwal Magadi

If and are functions of θ, and θ is the independent variable, then what exactly is Ψ? How does it look like on a graph? Can you sketch it?

The i equation is not an issue, but it depends on m and n, and they are merely binary signals, consisting of only 1's and 0's arranged in some order. Technically, they are just Heaviside step functions that switch when the conditions are met.

See example for m:

theta = linspace(-1, 1, 2001);
Psi   = theta;
Psi1  = sin(2*pi*Psi);
plot(theta, Psi, theta, Psi1), grid on
xlabel('\theta', 'fontsize', 16)
legend('\Psi', '\Psi_{1}', 'fontsize', 16, 'location', 'east')

m = 1/2*(sign(Psi - Psi1) + 1);     % the math of Heaviside step function
plot(theta, m, 'linewidth', 2), grid on, ylim([-0.5 1.5])
xlabel('\theta', 'fontsize', 16)
ylabel('m', 'fontsize', 16)

I think you get the idea for now. But the real Ψ does not necessarily behave like a straight line.

For the rest, you can follow @Jon's guidance.

[Prajwal Magadi]

Thank you very much for ur suggestion @Sam Chak.