Solve differential equation using Runge-Kutta

1 次查看(过去 30 天)
Hello everyone!
I have to solve the following equation by using the Runge-Kutta method:
all parameters are known, furthermore, Pin (function of t) and the timespan are vectors of dimensions 1x4096, so N(t) should be of the same dimension.
My code is the following:
n_dz=10;
dz=linspace(0,L,n_dz);
for i=2:n_dz
Dz=dz(i)-dz(i-1)
%%Compute N(t) using Runge-Kutte Method
Fun=@(N) (I/(q*V))-A*N-B*N.^2-C*N.^3-Gamma*(S*(N-N_tran)...
-k*(lambda-lambda_peak)^2).*(Pin/Dz)...
*(exp(((Gamma*(S*(N-N_tran)-k*(lambda-lambda_peak)^2))-alfa_s)*dz(i))...
-exp(((Gamma*(S*(N-N_tran)-k*(lambda-lambda_peak)^2))-alfa_s)*dz(i-1)))...
/((Gamma*(S*(N-N_tran)-k*(lambda-lambda_peak)^2))-alfa_s)...
*Dz/(V*h*f);
y=Runge_Kutta_4(time,Dt,N_start,Fun);
...
...
end
with Runge_Kutta_4 implemented as:
function y=Runge_Kutta_4(timespan,h,N_start,fun)
t=timespan;
n=length(t);
y=zeros(n,1);
y(1)=N_start;
for j=1:n-1
k1 = h*feval(fun,t(j),y(j));
k2 = h*feval(fun,t(j)+h/2,y(j)+k1/2);
k3 = h*feval(fun,t(j)+h/2,y(j)+k2/2);
k4 = h*feval(fun,t(j)+h,y(j)+k3);
y(j+1) = y(j)+1/6*(k1+2*k2+2*k3+k4);
end
end
Notice that in my code, the analysis over the total length "L" is divided into span of length "dz", therefore the last part of "Fun" is different from the formula shown above because it is due to an integral operation.
The following errors are shown:
??? Error using ==>
@(N)(I/(q*V))-A*N-B*N.^2-C*N.^3-Gamma*(S*(N-N_tran)-k*(lambda-lambda_peak)^2).*(Pin/Dz)*(exp(((Gamma*(S*(N-N_tran)-k*(lambda-lambda_peak)^2))-alfa_s)*dz(i))-exp(((Gamma*(S*(N-N_tran)-k*(lambda-lambda_peak)^2))-alfa_s)*dz(i-1)))/((Gamma*(S*(N-N_tran)-k*(lambda-lambda_peak)^2))-alfa_s)*Dz/(V*h*f)
Too many input arguments.
Error in ==> Runge_Kutta_4 at 7
k1 = h*feval(fun,t(j),y(j));
Error in ==> SOA_v1 at 47
y=Runge_Kutta_4(time,Dt,N_start,Fun);
Error in ==> RZ_source at 29
How can i solve this problem?
Thank you for your attention
Best regards

回答(1 个)

Tim Berk
Tim Berk 2017-9-19
I'm not sure I understand what your function (Fun) SHOULD be, but this is what you are telling it to be:
Fun = @(N) (I/(q*V))-A*N...
means Fun is a function of N, given by (I/(q*V))-A*N....
So you can put a value of N=9 into Fun using
feval(Fun,9)
And the output will be (I/(q*V))-A*9....
However, you are trying to put two values into the function by
k1 = h*feval(fun,t(j),y(j));
Hence the error Too many input arguments.

Community Treasure Hunt

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

Start Hunting!

Translated by