# Does anybody know how to calculate self-calling differential operators fast?

1 ビュー (過去 30 日間)
Mohammad Shojaei Arani 2021 年 11 月 10 日
コメント済み: Walter Roberson 2021 年 11 月 10 日
Hello friends,
I have a code with a for-loop where inside the for-loop there is a differential operator which uses a function, say A0, as input and performs
some operations and fthe results of the calculations are asigned to A0 again to be used for the next iteration. Intuitively, we understand that
the expressions in the for-loop get more and more complex at each iteration so that the computation time increases.
Long story short: In bellow I attach my code. The double for-loop (with tic and toc) becoms computationlly expensive
very when I increase J and K. I hope you have an idea to imrpve the speed of this code!
Babak
clc;
J=4;K=4;
H=sym('H',[J 1]);dt=0.1;
syms x par1 par2 par3 par4 par5
for j=1:J
H(j)=horner(simplify(exp(x^2/2)*diff(exp(-x^2/2),j)),x);
end
syms y y0 positive
muY(y)=-(y*(y*(y*(par1*y*par5^2 - par1*par2*par5) + par2*par3 + par1*par5^2) - par1*par2*par5))/(par2*par5*y^2 + par2*par5);
Hz=subs(H,x,(y-y0)/dt^(1/2));
eta=sym('eta',[J 1]);
tic;
for j=1:J
A0=Hz(j);S=A0;
for k=1:K
A0=muY(y)*diff(A0,y)+1/2*diff(A0,y,2);
S=S+A0*dt^k/factorial(k);
end
eta(j)=1/factorial(j)*limit(S,y,y0);
end
toc;
##### 1 件のコメント表示非表示 なし
Walter Roberson 2021 年 11 月 10 日
I wonder...
If you were to use
syms MUY(y)
and
A0=MUY(y)*diff(A0,y)+1/2*diff(A0,y,2);
then you should be able to get out a response in terms of powers of MUY(y) and derivatives of MUY(y) . At that point, you can use findSymType() to locate derivatives of MUY, and find the maximum derivative, as well as building up a list of which derivatives are used. Loop through iteratively building up the derivatives to the maximum used. Then subs() for the derivative actually used in S.
I think doing this should lower the cost of doing the derivatives of muY. The subs() at the end might be expensive, but I suspect you will still end up saving time.

サインインしてコメントする。

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by