Search All of the Math Forum:
Views expressed in these public forums are not endorsed by
Drexel University or The Math Forum.


Torsten
Posts:
1,516
Registered:
11/8/10


Re: Writing a system of differential equations correctly
Posted:
Apr 17, 2013 8:28 AM


"Douglas Nicolin" <douglas.nicolin@gmail.com> wrote in message <kkm2f3$2se$1@newscl01ah.mathworks.com>... > "Torsten" wrote in message <kklgfe$g4l$1@newscl01ah.mathworks.com>... > > > > > > I'm still dealing with errors when running the code I developed. Actually the error arise because of a mistake of writing in my opinion. But I'm able to see what's wrong until now. My function is shown below: > > > > > > function DXDt = testsol(t,X) > > > > > > global N > > > > > > N = 50; > > > r = 1; > > > dr = r/N; > > > Xeq = 1.7409; > > > rhoss = 1057; > > > rhoagua = 1000; > > > Do = 3e11; > > > k1 = 1e2; > > > > > > F(1) = 3*Do*exp(k1*X(1))/(X(N)^2)*((2*X(2)2*X(1))/(dr)^2); > > > > > > for i = 2:N2 > > > F(i) = Do*exp(k1*X(i))/(X(N)^2)*(((X(i+1)X(i1))/(2*dr))*... > > > ((2/(i*dr))+k1*(X(i+1)X(i1))/(2*dr))+(X(i+1)2*X(i)+X(i1))/(dr)^2)+... > > > (i*dr/X(N)^2)*(rhoss/rhoagua)*Do*exp(k1*X(N1))*((XeqX(N2))/dr)*((X(i+1)X(i1))/dr); > > > end > > > > > > for i = N1 > > > F(i) = Do*exp(abs(k1*X(i))/(X(N)))*(((XeqX(i1))/(2*dr))*... > > > ((2/(i*dr))+k1*(XeqX(i1))/(2*dr))+(Xeq2*X(i)+X(i1))/(dr)^2)+... > > > (i*dr/X(N)^2)*(rhoss/rhoagua)*Do*exp(k1*X(N1))*((XeqX(N2))/dr)^2; > > > end > > > > > > F(N) = ((rhoss/rhoagua)*Do*exp(k1*X(N1))/(X(N)))*(XeqX(N2))/dr; > > > > > > DXDt = F'; > > > > > > The main problem is that the last differential equation, defined in F(N) is not part of the system generated after using finite differences in the pde. This ode is the last one to be added that I mentioned before. However, the variable X(N) must represent the dependent variable of the last ode and it must appear in the system as well. When I run the whole code for the solution of my function I get this error message: > > > "??? Attempted to access X(50); index out of bounds because numel(X)=1." > > > which is natural because X(N) is not calculated before. My problem is that X(N) is a variable that must appear in the whole system, but is also a dependent variable to be calculated by integration in time. That's the biggest mistake I'm making: I'm not sure how to write this function properly. > > > > You must solve your problem in _one_ complete system of ODEs, not by calling the ODE integrator twice  once for X and once for R (that's what I think you are trying to do ?). > > Thus if the number of grid points is N, the dimension of the vector X and the dimension of the vector F from above must be N+1. Here, X(1) to X(N) represent your original solution variable X in the N grid points and element X(N+1) represents your solution variable R ; elements F(1) to F(N) represent the temporal derivatives of your solution variable X in the N grid points and element F(N+1) represents dR/dt. > > > > Best wishes > > Torsten. > > > You are wright. But in my system, the variable X(N) represents what would be my R(1) let's say. Once in the position N, in my case, there is the time derivative dRdt, I considered that the variable X(N) is my "R" variable that appears all over the system and in the last equation. That's my real problem. If my X(N) were a variable solve as a part of the system formed by my pde, there would be no errors on my code. But it's not. Adapting the indexes for the whole system is the problem. In other problems, when the dependence of the system on the X(N) variable does not exist it is easy to write the correct function. I'm search for examples similar to may case day by day and until now, I got nothing. I found many articles showing similar problems but solved in a different way. > > Kind regards > Douglas
Write the solution vector X in local variables:
function DXDt = testsol(t,X) global N N = 50; r = 1; dr = r/N; Xeq = 1.7409; rhoss = 1057; rhoagua = 1000; Do = 3e11; k1 = 1e2;
X_pde=zeros(N1);
For i=1:N1 X_pde(i)=X(i); end R=X(N); ...
Now you can work with X_pde and R. Does this solve your problem ?
Best wishes Torsten.



