Help me use this formula to complete this task (Runge–Kutta methods)

Question:

Use this method to solve the following IVP:

A series RC circuit with R = 5 W and C = 0.02 F is connected with a battery of E = 100 V. At t = 0, the voltage across the capacitor is zero.

(a) Obtain the subsequent voltage across the capacitor.

(b)find the true error if the exact solution is: Vc = V(1 – e^(-t/RC))

Hint:

Use the formula

r(dq/dt)+(q/c)=v

My work so far:

R=5;
C=0.02;
E=100;
Y0=0;
Xi=0;
Xf=1; 
h=1;
X=Xi:h:Xf;
Ydot=@(X,Y) Y*(1-exp(-X/R*C));




N=(Xf-Xi)/h;
for i=1:N
K1=Ydot(Xi,Y0);
K2=Ydot(Xi+0.5*h,Y0+0.5*h*K1);
K3=Ydot(Xi+ 1/2*h , Y0 +1/2*h*K2);
K4=Ydot(Xi+h, Y0+ h*K3);
Y(i+1)=Y0+h/6*(K1+2*K2+2*K3+K4);
Y0=Y(i+1);
end

I am a bit confused in the variables, how can I fix my code to work?

Answer 1

[James Tursa]

First, your derivative function should be obtained from your differential equation, not from the solution as you have it. You were given the following:

r(dq/dt)+(q/c)=v

So first solve this for the derivative (dq/dt):

(dq/dt) = (v - q/c)/r

Now you can use this expression to form a function handle for calculating this derivative given the state q:

v = something  <-- you fill this in
c = something  <-- you fill this in
r = something  <-- you fill this in
Ydot = @(t,q) (v-q/c)/r;  % based on the DE, not on the solution! Make sure you define v, c, r first!

(Side Note: I am using the (t,q) signature because it matches MATLAB syntax, even though your DE doesn't use t)

In your loop, you need to be calling the derivative function with the current state to calculate the derivative at that point. You are only doing this for Y currently, but this would only work for you in cases where the derivative doesn't depend on t. In your particular case this is true, but I would advise setting up your code so that it works regardless.E.g. something like this using X(i) on the rhs:

Y(1) = Y0;
for i=1:N
    K1 = Ydot(X(i)      ,Y(i)         );
    K2 = Ydot(X(i)+0.5*h,Y(i)+0.5*h*K1);
       :
    Y(i+1) = Y(i) + h/6*(K1+2*K2+2*K3+K4);
end

[Youssef Hani]

Hey James,

I am very thankful for your answer. It definitely cleared up some of the confusion I had. Below is my edited code according to the remarks you made.

Does this code correctly "Obtain the subsequent voltage across the capacitor" and "find the true error" for this particular case? Please let me know asap as I am expected to submit it tomorrow.

Again, much appreciated and thanks for your help.

R=5;
C=0.02;
V=100;
Y0=0;
Xi=0;
Xf=1; 
h=1;
X=Xi:h:Xf;
Ydot=@(t,q) (V-q/C)/R;
N=(Xf-Xi)/h;
for i=1:N
K1=Ydot(Xi,Y0);
K2=Ydot(Xi+0.5*h,Y0+0.5*h*K1);
K3=Ydot(Xi+ 1/2*h , Y0 +1/2*h*K2);
K4=Ydot(Xi+h, Y0+ h*K3);
Y(i+1)=Y0+h/6*(K1+2*K2+2*K3+K4);
Y0=Y(i+1);
end

Exact=V*(1-exp(-1/R*C));
Error=(Exact-Y0)/Exact;

[James Tursa]

You are very close. A couple of suggestions / corrections:

1) Drop the stepsize. You currently have h=1, but 1 is your final time so this gives you only one step. Make it quite a bit smaller, something like h = 0.01 would give you 100 steps.

h = 0.01;

2) Your exact formula is not quite right. Both the R and C need to be in the denominator. Also, you need to feed X into this exact formula (i.e., X is the same as t here). E.g.,

Vexact = V * (1 - exp(-X/(R*C)));

3) The q in the differential equation is the charge. So your Y variable is charge. But the question asks for voltage. To get the voltage you need to divide q (or Y in your case) by the capacitance. Then this is what gets compared to the exact voltage solution. So, after your for-loop is done,

Vrk4 = Y / C;

So now you can compare Vrk4 to Vexact. E.g., you could plot(X,Vrk4,X,Vexact), or plot(X,Vrk4-Vexact), etc.

4) If you want to get relative error and divide by the exact solution Vexact (to get a time history of the error), be sure to use the element-wise division operator ./ with the dot!

Again, I would remind you that your RK4 equations in the loop only work because your derivative function is not dependent on t. Your code would not work if the derivative function depended on t. So in subsequent assignments in your class, if you get a problem that has a t in the derivative function you will need to fix your code for this.

[Youssef Hani]

R=5;
C=0.02;
V=100;
Y0=0;
Xi=0;
Xf=1; 
h=0.01;
X=Xi:h:Xf;
Ydot=@(t,q) (V-q/C)/R;
N=(Xf-Xi)/h;
for i=1:N
K1=Ydot(Xi,Y0);
K2=Ydot(Xi+0.5*h,Y0+0.5*h*K1);
K3=Ydot(Xi+ 1/2*h , Y0 +1/2*h*K2);
K4=Ydot(Xi+h, Y0+ h*K3);
Y(i+1)=Y0+h/6*(K1+2*K2+2*K3+K4);
Y0=Y(i+1);
end
Vrk4= Y0 / C;
Exact=V * (1 - exp(-X/(R*C)));
Error=(Exact-Vrk4)./Exact;

In the 9th line "Ydot=@(t,q) (V-q/C)/R;" Is this correct? q is not previously declared.

[James Tursa]

"... q is not previously declared ..."

q and t are input arguments to this function. They do not have to be previously declared in your code.

Also, if you are only comparing the final voltage, then don't pass in the entire X to your Exact calculation as I stated. Instead, pass in only X(end).

Finally, you might make this adjustment for the N calculation to account for floating point effects:

N = round((Xf-Xi)/h);

To make your code more robust there are other things you could do to cover the cases where h doesn't divide into Xf-Xi evenly, but we will leave that for another time.

[Steven Lord]

One kind of stylistic suggestion for the future: the limit on the maximum length a variable name can have in MATLAB is currently 63 characters. Take advantage of this by giving your variables longer but more descriptive names. I'm not saying you need to use all 63 characters, but compare this:

Xi=0;
Xf=1; 
h=0.01;
X=Xi:h:Xf;

and this:

tinitial = 0;
tfinal = 1;
stepsize = 0.01;
tvector = tinitial:stepsize:tfinal;

I find the second version much easier to read and understand. To some extent, the second version documents itself as long as you know what the colon operator does.

Your assignment may mandate the name of the variable in which your script must store the final result. But that doesn't mean you can't write the code with descriptive names and define the required variable at the end. Consider what this line does if you've been storing your solution in a variable named solution.

Vrk4 = solution(end)/C;