Drexel dragonThe Math ForumDonate to the Math Forum



Search All of the Math Forum:

Views expressed in these public forums are not endorsed by Drexel University or The Math Forum.


Math Forum » Discussions » Software » comp.soft-sys.matlab

Topic: Writing a system of differential equations correctly
Replies: 6   Last Post: Apr 17, 2013 8:28 AM

Advanced Search

Back to Topic List Back to Topic List Jump to Tree View Jump to Tree View   Messages: [ Previous | Next ]
Torsten

Posts: 1,439
Registered: 11/8/10
Re: Writing a system of differential equations correctly
Posted: Apr 17, 2013 8:28 AM
  Click to see the message monospaced in plain text Plain Text   Click to reply to this topic Reply

"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 = 3e-11;
> > > k1 = 1e-2;
> > >
> > > F(1) = 3*Do*exp(k1*X(1))/(X(N)^2)*((2*X(2)-2*X(1))/(dr)^2);
> > >
> > > for i = 2:N-2
> > > F(i) = Do*exp(k1*X(i))/(X(N)^2)*(((X(i+1)-X(i-1))/(2*dr))*...
> > > ((2/(i*dr))+k1*(X(i+1)-X(i-1))/(2*dr))+(X(i+1)-2*X(i)+X(i-1))/(dr)^2)+...
> > > (i*dr/X(N)^2)*(rhoss/rhoagua)*Do*exp(k1*X(N-1))*((Xeq-X(N-2))/dr)*((X(i+1)-X(i-1))/dr);
> > > end
> > >
> > > for i = N-1
> > > F(i) = Do*exp(abs(k1*X(i))/(X(N)))*(((Xeq-X(i-1))/(2*dr))*...
> > > ((2/(i*dr))+k1*(Xeq-X(i-1))/(2*dr))+(Xeq-2*X(i)+X(i-1))/(dr)^2)+...
> > > (i*dr/X(N)^2)*(rhoss/rhoagua)*Do*exp(k1*X(N-1))*((Xeq-X(N-2))/dr)^2;
> > > end
> > >
> > > F(N) = ((rhoss/rhoagua)*Do*exp(k1*X(N-1))/(X(N)))*(Xeq-X(N-2))/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 = 3e-11;
k1 = 1e-2;

X_pde=zeros(N-1);

For i=1:N-1
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.




Point your RSS reader here for a feed of the latest messages in this topic.

[Privacy Policy] [Terms of Use]

© Drexel University 1994-2014. All Rights Reserved.
The Math Forum is a research and educational enterprise of the Drexel University School of Education.