How to Implement 3rd order equation in matlab script?
4 次查看(过去 30 天)
显示 更早的评论
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.
1 个评论
Sam Chak
2023-7-25
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)
采纳的回答
Jon
2023-7-24
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)
6 个评论
Jon
2023-7-25
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
2023-7-27
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?
更多回答(1 个)
Sam Chak
2023-7-25
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.
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Calculus 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!