- Do you receive warning and/or error messages? If so the full and exact text of those messages (all the text displayed in orange and/or red in the Command Window) may be useful in determining what's going on and how to avoid the warning and/or error.
- Does it do something different than what you expected? If so, what did it do and what did you expect it to do?
- Did MATLAB crash? If so please send the crash log file (with a description of what you were running or doing in MATLAB when the crash occured) to Technical Support so we can investigate.
이 질문을 팔로우합니다.
- 팔로우하는 게시물 피드에서 업데이트를 확인할 수 있습니다.
- 정보 수신 기본 설정에 따라 이메일을 받을 수 있습니다.
I want to use fzero to find the root of a transcendental function r. I have problem on line 33-37 and 52-63
조회 수: 1 (최근 30일)
이전 댓글 표시
%% initialization
clear
clc
close all
%% Model parameters
global k1 eta1 alpha3 gamma1 d N Phi xi h C q qstar Theta OddAsymptotes m ...
qq r1 r2 alpha aa b Efun F B Ainv r v theta_sol Vfun v_sol p ModifiedOddAsym ...
qstarTemporary theta_init z
k1 = 6*10^(-12); % elastic constant
eta1 = 0.0240; % viscosity
xi = 0.585; % activity parameter
alpha3 = -0.001104; % viscosity
gamma1 = 0.1093; % viscosity
Theta = 0.0001;
theta_init = Theta.*sin(pi*z/d);
d = 0.0002;
N = 20;
% Simplified parameters
alpha=1-alpha3^2/gamma1/eta1;
C = (alpha3*xi*alpha/eta1)./(4*k1*q.^2/d^2-alpha3*xi/eta1);
% Plotting r(q)
q=0:0.01:(5*pi);
r=q-(1-alpha).*tan(q)+ C.*tan(q);
Error using .*
Arrays have incompatible sizes for this operation.
Arrays have incompatible sizes for this operation.
figure
plot(q,r)
axis([0 5*pi -20 20])
% Assign g for q and plot r, and compute the special asymptotes
syms C
for q=0:N
qstar(i) = solve(1/C == 0);
end
% Compute Odd asymptote
for j=1:N
OddAsymptotes(j)= 1/2*(2*j + 1)*pi;
end
ModifiedOddAsym = [abs(OddAsymptotes) abs(qstar)];
qstarTemporary = min(ModifiedOddAsym);
% Compute the rest of the Roots
for j=1:N
OddAsymptotes(j)= 1/2.*(2*j + 1)*pi;
end
% Compute all symptotes
AllAsymptotes = [sort(OddAsymptotes) sort(qstar)];
q0 =0;
for q = 0.001:AllAsymptotes(1);
r1 = fzero(r, q0)
end
for i=1:N
for q =AllAsymptotes(i):AllAsymptotes(i+1);
r2 = fzero(r, q0)
end
end
% Compute all roots
gg = [r1 r2];
% Remove the repeated roots if any & redefine n
qq = unique(gg);
m = length(qq);
% Compute the time constants tan_n
for i=1:m
p(i) = (gamma1*alpha)./(4*k1*qq(i)^2/d^2-alpha3*xi/eta1);
end
% Compute integral coefficients for off-diagonal elements θ_n matrix
for i=1:m
for j= i+1:m
for z=0:d
b(i, j) = int(Efun(i).*Efun(j));
b(j, i) = b(i, j);
aa(i, j) = b(i, j);
aa(j, i)= b(j, i);
end
end
end
% Calculate integral coefficients for diagonal elements in theta_n matrix
for i=1:m
for z=0:d
aa(i, i)= int(Efun(i).^2);
end
end
% Calculate integrals for RHS vector
for i=1:m
for z=0:d
F(i) = int(theta_init.*Efun(i));
end
end
% Create matrix A and RHS vector B
for i=1:m
for j= 1:m
aa(i, j)
end
end
for i=1:m
B = F(i)
end
% Calculate inverse of A and solve A*x=B
Ainv = inv(A);
r = B*Ainv;
% Define Theta(z,t)
for i=1:m
theta_sol(i) = sum(r(i).*Efun(i).*exp(-t/p(i)));
end
% Compute v_n for n times constant
for i=1:m
v(i) = (-2*k1*alpha3*qq(i))*(1/(d^2*eta1*alpha*gamma1))+ alpha3^2*xi/(2*(eta1)^2*qq(i)*alpha*gamma1)+xi/(2*eta1*qq(i));
end
% Compute v(z,t) from initial conditions
for i=1:m
Vfun(i) = d*sin(qq(i)-2*qq(i)*z/d)+(2*z-d)*sin(qq(i));
end
% Define v(z,t)
for i=1:m
v_sol = sum(v(i).*r(i)*Vfun(i).*exp(-t/p(i)));
end
댓글 수: 17
Steven Lord
2022년 8월 23일
What does "have problem" mean in this context?
Also please indicate specifically which lines are lines 33-37 and 52-63 in the code you posted above.
University Glasgow
2022년 8월 23일
I receive the following erros:
Array indices must be positive integers or logical values. Error in Case1_Ijuptilk (line 36)
qstar(i) = solve(1/C == 0);
Error using fzero :FUN must be a function, a valid character vector expression, or an inline function object. Error in Case1_Ijuptilk (line 56)
r1 = fzero(r, q0);
Error using fzero :FUN must be a function, a valid character vector expression, or an inline function object. Error in Case1_Ijuptilk (line 61)
r1 = fzero(r, q0);
Torsten
2022년 8월 23일
Array indices must be positive integers or logical values. Error in Case1_Ijuptilk (line 36)
qstar(i) = solve(1/C == 0);
i must be q. But what variable do you want to solve for ? C is quadratic in q - so it permits at most two solutions.
r1 = fzero(r, q0);
r in your code is an array of values for q =0:0.01:(5*pi). "fzero" needs a function as first argument.
Define
alpha=1-alpha3^2/gamma1/eta1;
C = @(q) (alpha3*xi*alpha/eta1)./(4*k1*q.^2/d^2-alpha3*xi/eta1);
r = @(q) q-(1-alpha).*tan(q)+ C(q).*tan(q);
University Glasgow
2022년 8월 23일
편집: University Glasgow
2022년 8월 23일
I want to solve for 1/C for q from 0 to N. I tried something like ths
C = @(q) (alpha3*xi*alpha/eta1)./(4*k1*q.^2/d^2-alpha3*xi/eta1);
for q=1:N
qstar(q) = solve(1/C);
end
Torsten
2022년 8월 23일
I don't understand what you mean by "I want to solve for 1/C". Do you mean: "I want to evaluate 1/C for q = 1,...,N" ?
University Glasgow
2022년 8월 23일
Thank you for your help. Please hold on for now. I'm new to matlab.
Can help me to check why the code below print same values of q
d = 0.2;
N =10;
fun = @(x, d) d.*sin((i)*x);
for i=1:N
q = integral(@(x) fun(x, d) ,0,1,'ArrayValued',true)
end
University Glasgow
2022년 8월 23일
I have revised the code but I'm still receiving an error: Arrays have incompatible sizes for this operation
Error in Case1_Ijuptilk (line 28) r=q-(1-alpha).*tan(q)+ C.*tan(q);
%% initialization
clear
clc
close all
%% Model parameters
global k1 eta1 alpha3 gamma1 d N xi C q qstar Theta OddAsymptotes m ...
qq r1 r2 alpha aa b Efun F B Ainv r v theta_sol Vfun v_sol p ModifiedOddAsym ...
qstarTemporary theta_init z AllAsymptotes gg
k1 = 6*10^(-12); % elastic constant
eta1 = 0.0240; % viscosity
xi = 0.01; % activity parameter
alpha3 = -0.001104; % viscosity
gamma1 = 0.1093; % viscosity
Theta = 0.0001;
theta_init = Theta.*sin(pi*z/d);
d = 0.0002;
N = 20;
% Simplified parameters
alpha=1-alpha3^2/gamma1/eta1;
C = (alpha3*xi*alpha/eta1)./(4*k1*q.^2/d^2-alpha3*xi/eta1);
% Plotting r(q)
q=0:0.01:(5*pi);
r=q-(1-alpha).*tan(q)+ C.*tan(q);
figure
plot(q,r)
axis([0 5*pi -20 20])
% Assign g for q and plot r, and compute the special asymptote
for q=1:N
qstar(q) = 1./(alpha3*xi*alpha/eta1)./(4*k1*q.^2/d^2-alpha3*xi/eta1);
end
% Compute Odd asymptote
for j=1:N
OddAsymptotes(j)= 1/2*(2*j + 1)*pi;
end
ModifiedOddAsym = [abs(OddAsymptotes) abs(qstar)];
qstarTemporary = min(ModifiedOddAsym);
% Compute the rest of the Roots
for j=1:N
OddAsymptotes(j)= 1/2.*(2*j + 1)*pi;
end
% Compute all symptotes
r = @(q, C) q-(1-alpha).*tan(q)+ C.*tan(q);
AllAsymptotes = [sort(OddAsymptotes) sort(qstar)];
q0 =0;
for q = 1:AllAsymptotes(1)
r1(q) = fzero(@(q)r(q,C), q0);
end
for i=1:N
for q = 2:N
r2(q) = fzero(@(q)r(q,C), q0);
end
end
% Compute all roots
gg = [r1 r2];
% Remove the repeated roots if any & redefine n
qq = unique(gg);
m = length(qq);
% Compute the time constants tan_n
for i=1:m
p(i) = (gamma1*alpha)./(4*k1*qq(i)^2/d^2-alpha3*xi/eta1);
end
% Compute theta_n from initial conditions
Efun = @(z, i, d)cos(qq(i)-2*qq(i).*z/d)-cos(qq(i));
% Compute integral coefficients for off-diagonal elements θ_n matrix
for i=1:m
for j= i+1:m
b(i, j) = integral(@(x)Efun(x, i, d).*@(x)Efun(x, j, d), 0, d);
b(j, i) = b(i, j);
aa(i, j)= b(i, j);
aa(j, i)= b(j, i);
end
end
% Calculate integral coefficients for diagonal elements in theta_n matrix
for i=1:m
aa(i, i)= integral((@(x)Efun(x, i, d)).^2, 0, d);
end
% Calculate integrals for RHS vector
for i=1:m
F(i) = integral(theta_init.*@(x)Efun(i, x, d), 0, d);
end
% Create matrix A and RHS vector B
for i=1:m
for j= 1:m
aa(i, j);
end
end
for i=1:m
B = F(i);
end
% Calculate inverse of A and solve A*x=B
Ainv = inv(A);
r = B.*Ainv;
% Define Theta(z,t)
for i=1:m
theta_sol(i) = sum(r(i).*@(x)Efun(i, x, d).*exp(-t/p(i)));
end
% Compute v_n for n times constant
for i=1:m
v(i) = (-2*k1*alpha3*qq(i))*(1/(d^2*eta1*alpha*gamma1))+ alpha3^2*xi/(2*(eta1)^2*qq(i)*alpha*gamma1)+xi/(2*eta1*qq(i));
end
% Compute v(z,t) from initial conditions
for i=1:m
Vfun(i) = d*sin(qq(i)-2*qq(i)*z/d)+(2*z-d)*sin(qq(i));
end
% Define v(z,t)
for i=1:m
v_sol = sum(v(i).*r(i)*Vfun(i).*exp(-t/p(i)));
end
University Glasgow
2022년 8월 23일
편집: Torsten
2022년 8월 24일
Thank you for your support. This revised version. I want to use sum function to compute for theta_sol on line 113 and vsol on line 125. I expect to get vectors then plot them but theta_sol and vsol is 0.
%% initialization
clear
clc
close all
%% Model parameters
global k1 eta1 alpha3 gamma1 d N xi C q qstar Theta OddAsymptotes m ...
qq alpha aa b Efun F B r p ModifiedOddAsym qstarTemporary ...
AllAsymptotes gg Efun2 Efun1 v theta_sol Vfun Efun3 v_sol
k1 = 6*10^(-12); % elastic constant
eta1 = 0.0240; % viscosity
xi = 0.01; % activity parameter
alpha3 = -0.001104; % viscosity
gamma1 = 0.1093; % viscosity
Theta = 0.0001;
d = 0.0002;
N = 20;
% Simplified parameters
alpha=1-alpha3^2/gamma1/eta1;
C = (alpha3*xi*alpha/eta1)./(4*k1*q.^2/d^2-alpha3*xi/eta1);
Unrecognized function or variable 'q'.
%% Plotting r(q)
t = [0:0.1:10];
nt=length(t);
q=0:0.01:(5*pi);
r=q-(1-alpha).*tan(q)+ C.*tan(q);
figure
plot(q,r)
axis([0 5*pi -20 20])
% Assign g for q and plot r, and compute the special asymptote
r=q-(1-alpha).*tan(q)+ C.*tan(q);
for q=1:N
qstar(q) = 1./(alpha3*xi*alpha/eta1)./(4*k1*q.^2/d^2-alpha3*xi/eta1);
end
% Compute Odd asymptote
for j=1:N
OddAsymptotes(j)= 1/2*(2*j + 1)*pi;
end
ModifiedOddAsym = [abs(OddAsymptotes) abs(qstar)];
qstarTemporary = min(ModifiedOddAsym);
% Compute the rest of the Roots
for j=1:N
OddAsymptotes(j)= 1/2.*(2*j + 1)*pi;
end
% Compute all symptotes
AllAsymptotes = [sort(OddAsymptotes) sort(qstar)];
gg = AllAsymptotes;
% Remove the repeated roots if any & redefine n
qq = unique(gg);
m = length(qq);
% Compute the time constants tan_n
for i=1:m
p(i) = (gamma1*alpha)./(4*k1*qq(i)^2/d^2-alpha3*xi/eta1);
end
% Compute theta_n from initial conditions
Efun = @(z, i, d)cos(qq(i)-2*qq(i).*z/d)-cos(qq(i));
Efun2 = @(z, j, d)cos(qq(i)-2*qq(j).*z/d)-cos(qq(j));
% Compute integral coefficients for off-diagonal elements θ_n matrix
b = zeros(m, m);
for i=1:m
for j= i+1:m
b(i, j) = integral(@(z)(Efun(z, i, d).*Efun2(z, j, d)), 0, d);
b(j, i) = b(i, j);
aa(i, j) = b(i, j);
aa(j, i) = b(j, i);
end
end
% Calculate integral coefficients for diagonal elements in theta_n matrix
for i=1:m
aa(i, i)= integral((@(z)Efun(z, i, d).^2), 0, d);
end
% Calculate integrals for RHS vector
Efun1 = @(z, i, d, Theta) Theta.*sin(pi*z/d).*cos(qq(i)-2*qq(i).*z/d)-cos(qq(i));
for i=1:m
F(i) = integral(@(z)Efun1(z, i, d, Theta), 0, d);
end
% Create matrix A and RHS vector BN_mat = N-2;
for i=1:m
for j= 1:m
A(i, j) = aa(i, j);
B(i) = F(i);
end
end
% Calculate inverse of A and solve A*x=B
Ainv = inv(A);
r = B.*Ainv;
% Compute theta(z,t) from initial conditionsDefine Theta(z,t)
for i =1:m
for z=0:d
Efun3(i) = cos(qq(i)-2*qq(i).*z/d)-cos(qq(i));
end
end
for i=1:m
theta_sol(i) = sum(r(i).*Efun3(i).*exp(-t/p(i)));
end
% Compute v_n for n times constant
for i=1:m
v(i) = (-2*k1*alpha3*qq(i))*(1/(d^2*eta1*alpha*gamma1))+ alpha3^2*xi/(2*(eta1)^2*qq(i)*alpha*gamma1)+xi/(2*eta1*qq(i));
end
% Compute v(z,t) from initial conditions and define v(z,t)
for i=1:m
for z=0:d
Vfun(i) = d*sin(qq(i)-2*qq(i)*z/d)+(2*z-d)*sin(qq(i));
v_sol = sum(v(i).*r(i)*Vfun(i).*exp(-t/p(i)));
end
end
University Glasgow
2022년 8월 25일
But i want to solve for the theta_v and v_sol for various times for z=0:d. I can help out please.
Walter Roberson
2022년 8월 25일
Remember that the defined increment is 1, so z=0:d is:
- empty if d is empty
- empty if d(1)<0
- [0] if d(1) < 1
- [0 1] if 1 <= d(1) < 2
- [0 1 2] if 2 <= d(1) < 3
and so on, integers only.
WIth your d = 0.0002, the d(1)<1 case applies, and you would get only [0]
University Glasgow
2022년 8월 25일
편집: Walter Roberson
2022년 8월 26일
This is the code I tried to translate to matlab. This is actually maple code but maple is limited at this point. I tried to use codegenerator in maple to translate the code but it didn't work. I'm entirely new to matlab. I started using matlab two weeks ago. Please can someone help me to fix the code using my code above.
doCalc:= proc( xi )
# Import Packages
uses ArrayTools, Student:-Calculus1, LinearAlgebra,
ListTools, RootFinding, plots, ListTools:
local gamma__1:= .1093,
alpha__3:= -0.1104e-2,
k__1:= 6*10^(-12),
d:= 0.2e-3,
theta0:= 0.001,
eta__1:= 0.240e-1,
alpha:= 1-alpha__3^2/(gamma__1*eta__1),
c:= alpha__3*xi*alpha/(eta__1*(4*k__1*q^2/d^2-alpha__3*xi/eta__1)),
theta_init:= theta0*sin(Pi*z/d),
n:= 30,
g, f, b1, b2, qstar, OddAsymptotes, ModifiedOddAsym,
qstarTemporary, indexOfqstar2, qstar2, AreThereComplexRoots,
soln1, soln2, qcomplex1, qcomplex2, gg, qq, m, pp, j, i,
AllAsymptotes, p, Efun, b, aa, F, A, B, Ainv, r, theta_sol, v, Vfun, v_sol,minp,nstar;
# Assign g for q and plot g, Set q as a complex and Compute the Special Asymptotes
g:= q-(1-alpha)*tan(q)+ c*tan(q):
f:= subs(q = x+I*y, g):
b1:= evalc(Re(f)) = 0:
b2:= evalc(Im(f)) = 0:
qstar:= (fsolve(1/c = 0, q = 0 .. infinity)):
OddAsymptotes:= Vector[row]([seq(evalf(1/2*(2*j + 1)*Pi), j = 0 .. n)]);
# Compute Odd asymptote
ModifiedOddAsym:= abs(`-`~(OddAsymptotes, qstar));
qstarTemporary:= min(ModifiedOddAsym);
indexOfqstar2:= SearchAll(qstarTemporary, ModifiedOddAsym);
qstar2:= OddAsymptotes(indexOfqstar2);
# Compute complex roots
AreThereComplexRoots:= type(true, 'truefalse');
try
#print({min(qstar2, qstar), max(qstar2, qstar)});
#print({eval(subs({x=min(qstar2, qstar)+0.1,y=0},b1))});
soln1:= fsolve({b1, b2}, {x = min(qstar2, qstar)+0.1 .. max(qstar2, qstar)-0.1, y = 0.1 .. 100});
soln2:= fsolve({b1, b2}, {x = min(qstar2, qstar) .. max(qstar2, qstar), y = -infinity .. 0});
#print(soln1);
qcomplex1:= subs(soln1, x+I*y);
qcomplex2:= subs(soln2, x+I*y);
catch:
AreThereComplexRoots:= type(FAIL, 'truefalse');
end try;
# Compute the rest of the Roots
OddAsymptotes:= Vector[row]([seq(evalf((1/2)*(2*j+1)*Pi), j = 0 .. n)]);
AllAsymptotes:= sort(Vector[row]([OddAsymptotes, qstar]));
if AreThereComplexRoots
then gg:= [ qcomplex1, qcomplex2,op(Roots(g, q = 0.1e-3 .. AllAsymptotes[1])),
seq(op(Roots(g, q = AllAsymptotes[i] .. AllAsymptotes[i+1])), i = 1 .. n)];
elif not AreThereComplexRoots
then gg:= [op(Roots(g, q = 0.1e-3 .. AllAsymptotes[1])), seq(op(Roots(g, q = AllAsymptotes[i] .. AllAsymptotes[i+1])), i = 1 .. n)];
end if:
# Remove the repeated roots if any & Redefine n
qq:= MakeUnique(gg):
m:= numelems(qq):
## Compute the
`τ_n`;
time constants
for i to m do
p[i] := gamma__1*alpha/(4*k__1*qq[i]^2/d^2-alpha__3*xi/eta__1):
end do:
## Compute θ_n from initial conditions
for i to m do
Efun[i] := cos(qq[i]-2*qq[i]*z/d)-cos(qq[i]):
end do:
## Compute integral coefficients for off-diagonal elements θ_n matrix
printlevel := 2:
for i to m do
for j from i+1 to m do
b[i, j] := int(Efun[i]*Efun[j], z = 0 .. d):
b[j, i] := b[i, j]:
aa[i, j] := b[i, j]:
aa[j, i] := b[j, i]:
end do :
end do:
## Calculate integral coefficients for diagonal elements in theta_n matrix
for i to m do
aa[i, i] := int(Efun[i]^2, z = 0 .. d):
end do:
## Calculate integrals for RHS vector
for i to m do
F[i] := int(theta_init*Efun[i], z = 0 .. d):
end do:
## Create matrix A and RHS vector B
A := Matrix([seq([seq(aa[i, j], i = 1 .. m)], j = 1 .. m)]):
B := Vector([seq(F[i], i = 1 .. m)]):
## Calculate inverse of A and solve A*x=B
Ainv := 1/A:
r := MatrixVectorMultiply(Ainv, B):
## Define Theta(z,t)
theta_sol := add(r[i]*Efun[i]*exp(-t/p[i]), i = 1 .. m):
## Compute v_n for n times constant
for i to m do
v[i] := (-2*k__1*alpha__3*qq[i])*(1/(d^2*eta__1*alpha*gamma__1))+ alpha__3^2*xi/(2*(eta__1)^2*qq[i]*alpha*gamma__1)+xi/(2*eta__1*qq[i]):
end do:
## Compute v(z,t) from initial conditions
for i to m do
Vfun[i] := d*sin(qq[i]-2*qq[i]*z/d)+(2*z-d)*sin(qq[i]):
end do:
## Define v(z,t)
v_sol := add(v[i]*r[i]*Vfun[i]*exp(-t/p[i]), i = 1 .. m):
##
minp:=min( seq( Re(p[i]), i=1..m) ):
member(min(seq( Re(p[i]), i=1..m)),[seq( Re(p[i]), i=1..m)],'nstar'):
## Return all the plots
return theta_init, theta_sol, v_sol, minp, eval(p), nstar, p[nstar], g, qq, qstar, AreThereComplexRoots;
end proc:
Steven Lord
2022년 8월 25일
You've posted a lot of code. But let's take a step back for a second. Can you explain in words not code the problem you're trying to solve? What do you know and what are you trying to find or compute? Knowing the context may help us offer better guidance in how to adjust the code.
University Glasgow
2022년 8월 26일
편집: Walter Roberson
2022년 8월 26일
First of all, I have a problem with computing the root of the function r. It seems that maple do skip some roots of the function r that's why I have change to matlab. I want to compute the roots of r for q=(2*i-2)*Pi/2..(2*i+1)*Pi/2) and assign the results to gg, if there is any repeated, i used unique() to remove the repeated roots and finally assign the result to qq and use qq to compute p(i), theta_sol and v_sol. I don't know how to find roots with matlab. I tried vpasolve() but i need syms which didn't work in my case.
Finally, I want to plot theta_sol and v_sol for various times at z=0:d.
%% initialization
clear
clc
close all
%% Model parameters
global k1 eta1 alpha3 gamma1 N xi C q qstar Theta OddAsymptotes m ...
qq alpha aa b Efun F B r p ModifiedOddAsym qstarTemporary ...
AllAsymptotes gg Efun2 Efun1 v theta_sol Vfun Efun3 v_sol d z
k1 = 6*10^(-12); % elastic constant
eta1 = 0.0240; % viscosity
xi = 0.01; % activity parameter
alpha3 = -0.001104; % viscosity
gamma1 = 0.1093; % viscosity
Theta = 0.0001;
d = 0.0002;
N = 20;
% Simplified parameters
alpha=1-alpha3^2/gamma1/eta1;
C = (alpha3*xi*alpha/eta1)./(4*k1*q.^2/d^2-alpha3*xi/eta1);
%% Plotting r(q)
t = [0:0.1:10];
nt=length(t);
q=0:0.01:(5*pi);
r=q-(1-alpha).*tan(q)+ C.*tan(q);
figure
plot(q,r)
axis([0 5*pi -20 20])
% Remove the repeated roots if any & redefine n
qq = unique(gg);
m = length(qq);
% Compute the time constants tan_n
for i=1:m
p(i) = (gamma1*alpha)./(4*k1*qq(i)^2/d^2-alpha3*xi/eta1);
end
% Compute theta_n from initial conditions
Efun = @(z, i, d)cos(qq(i)-2*qq(i).*z/d)-cos(qq(i));
Efun2 = @(z, j, d)cos(qq(i)-2*qq(j).*z/d)-cos(qq(j));
% Compute integral coefficients for off-diagonal elements θ_n matrix
b = zeros(m, m);
for i=1:m
for j= i+1:m
b(i, j) = integral(@(z)(Efun(z, i, d).*Efun2(z, j, d)), 0, d);
b(j, i) = b(i, j);
aa(i, j) = b(i, j);
aa(j, i) = b(j, i);
end
end
% Calculate integral coefficients for diagonal elements in theta_n matrix
for i=1:m
aa(i, i)= integral((@(z)Efun(z, i, d).^2), 0, d);
end
% Calculate integrals for RHS vector
Efun1 = @(z, i, d, Theta) Theta.*sin(pi*z/d).*cos(qq(i)-2*qq(i).*z/d)-cos(qq(i));
for i=1:m
F(i) = integral(@(z)Efun1(z, i, d, Theta), 0, d);
end
% Create matrix A and RHS vector BN_mat = N-2;
for i=1:m
for j= 1:m
A(i, j) = aa(i, j);
B(i) = F(i);
end
end
% Calculate inverse of A and solve A*x=B
Ainv = inv(A);
b = Ainv\B';
% Compute theta(z,t) from initial conditionsDefine Theta(z,t)
for i=1:m
z = 0:d;
Efun3(i) = cos(qq(i)-2*qq(i).*z/d)-cos(qq(i));
theta_sol(i) = sum(b(i).*Efun3(i).*exp(-t/p(i)));
end
% Compute v_n for n times constant
for i=1:m
v(i) = (-2*k1*alpha3*qq(i))*(1/(d^2*eta1*alpha*gamma1))+ alpha3^2*xi/(2*(eta1)^2*qq(i)*alpha*gamma1)+xi/(2*eta1*qq(i));
end
% Compute v(z,t) from initial conditions and define v(z,t)
for i=1:m
Vfun(i) = d*sin(qq(i)-2*qq(i)*z/d)+(2*z-d)*sin(qq(i));
v_sol = sum(v(i).*r(i)*Vfun(i).*exp(-t/p(i)));
end
답변 (0개)
참고 항목
카테고리
Help Center 및 File Exchange에서 Function Handles에 대해 자세히 알아보기
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!오류 발생
페이지가 변경되었기 때문에 동작을 완료할 수 없습니다. 업데이트된 상태를 보려면 페이지를 다시 불러오십시오.
웹사이트 선택
번역된 콘텐츠를 보고 지역별 이벤트와 혜택을 살펴보려면 웹사이트를 선택하십시오. 현재 계신 지역에 따라 다음 웹사이트를 권장합니다:
또한 다음 목록에서 웹사이트를 선택하실 수도 있습니다.
사이트 성능 최적화 방법
최고의 사이트 성능을 위해 중국 사이트(중국어 또는 영어)를 선택하십시오. 현재 계신 지역에서는 다른 국가의 MathWorks 사이트 방문이 최적화되지 않았습니다.
미주
- América Latina (Español)
- Canada (English)
- United States (English)
유럽
- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)
- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom (English)
아시아 태평양
- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)