Values in code are not updating ? /Not getting desired output.

I am trying to solve some couple ODE`s using runga-kutta 4th order NUmerical method,i am not getting desired result.

Answer 1

[Walter Roberson]

How can you tell you are not getting the desired result? You are not saving the results of your calculation, except that your array named "temporary" will reflect the result of the last for loop iteration (as you overwrite it on every iteration.)

[KAPIL MAMTANI]

Hi,Mr.Walter Thank you for responding. I previously tried saving results,for example in case of N in for as: N(i+1) = N(i) + (h/6)*(x1+2*x2+2*x3+x4); but it showed error ' In an assignment A(I) = B, the number of elements in B and I must be the same.', that is why i switched to using temporary as storage.Could you please tell me how to act now.

[Walter Roberson]

You have the code

 Gp = ((Kt*(1.015^(T-20)))*(P/(P2+P))*(TIC/(C2+TIC))*(N/N2+N));

Your vector P has a / operation involving something the same length as P so you are doing a least square fit that will end up with a scalar. I doubt that is what you want to do but at least it gives you a scalar. Likewise TIC/ something the same length as TIC will give you a scalar. But then you have your (N/2+N) where N is a vector. This is the same as N*1.5 and so gives a vector. Therefore your Gp gives a vector. You use Gp in constructing many of your other variables so most of your variables end up as vectors. That is why you cannot assign that result to N(i+1), because that is a scalar and you are trying to store a vector in it.

You need to learn how to use the debugger. Step through your code and see what the size of each result is and see if it makes sense to you.

[Walter Roberson]

Please notice that your X, Y, Z, and W definitions all ignore the t values being passed in. Is there a reason to pass the t values in at all? All of your calculations ignore the value of t, caring only how long it is.

[KAPIL MAMTANI]

First of all, 1.All the four variables i.e., TIC N A P are all having same definitions/type/class/size in workspace, so i could not understand why N/value+N is giving error and not other three. Also is there a way to convert N/N2+N into scalar before using it into Gp. 2. In x y z w ,the derivative is in respect of t as in d(N)/dt = RHS as defined, so t is utmost important to be passed.

[KAPIL MAMTANI]

I tried something like converting Gp to scalar value.the program is running and updating values but Gp value while being converted into scalar is using single value,and it is not being updating.I want Gp to be updated in each loop need assistance

[Walter Roberson]

You had expressions of the form

A/(A+B)

for vector A. The / operator is the mrdivide operator, a matrix algebra operation. Just the same way that A*B is defined as matrix multiplication in MATLAB, / is matrix division. When you have matrix division of two row vectors the same size the result is a scalar.

However you have N/N2 where N2 is scalar. This is a different case than the ones where the top and bottom were both vectors. You should consider the possibility that you intended to write N/(N2+N)

[Walter Roberson]

Look at your W definition

W = @(t,A) ((Ra-(Ka*neta)*A)+(CO2sat));

If it is of utmost importance that t be passed in then why is there no t in what is calculated? And is it t that is needed or is it dt that is needed?

[Walter Roberson]

If Gp should be updated in each loop then what should it be updated to and how should that affect the calculation? You use the norm of Gp now, so are you wanting the norm of the updated Gp?

It is not clear to me that you would want to use all of P and similar variables when you calculate Gp. It would make more sense to me if you were using only the values for the current time, or else only the values for the times from the beginning until the present.

[KAPIL MAMTANI]

Thanks for replying to me. Yeah that is exactly what i intended,for N.About W what i intended here is that RHS of equation d(TIC)/dt = ((Ra-(Ka*neta)*A)+(CO2sat)) , so i guess that is why i had to pass t here as TIC is a function of t and A. what say. I am attaching latest modified code file.

[Walter Roberson]

None of your formula appear to involve time explicitly. Instead you have the implicit that if d(TIC)/dt = ((Ra-(Ka*neta)*A)+(CO2sat)) then TIC = integral (Ra-(Ka*neta)*A)+(CO2sat)) with respect to time, and you are doing a particular form of numeric integration. Because the formula do not involve time specifically the numeric integrations turn out to only need dt, which is the same as your h, which you do use in the loop. But that doesn't mean that your W, X, Y, Z formula need to know t.

I am becoming more and more certain that your Gp should be a function rather than a constant.

[KAPIL MAMTANI]

• Thanks again.* I am also thinking of making Gp a function rather than a constant.So,you are saying that I can opt out t in w x y z definitions....

[Walter Roberson]

Right, you can remove t from your W X Y Z definitions, and do not need to pass the corresponding parameter when calling those functions.

[KAPIL MAMTANI]

Thank you very much Mr walter for all your help and guidance.

[KAPIL MAMTANI]

Hi, I have been running my code with different initial values for many runs but it is very rigorous work to find the optimum set of values.I want to check my code for varying input parameters/- Q(varying monthly),N(1),P(1),TIC(1). Also, I want to put a conditional statement for N/P and C/N in the code to change the formula as per scenario. Can you guide me in right direction?