My Numerical Solution doesn't align with the exact solution
显示 更早的评论
This is a 2nd order differential equation and I'm using ODE45 for my numerical Solution. My Numerical Solution is way too high. Notice that the values for numerical is *10^8. Additionally, I first solve my constants in order to make my equation code at matlab simpler.
My initial conditions are vc(0) = 11; vc' = 77459.66692
Exact Solution:
clear all
close all
clc
C = 3*10^(-6);
R = 12;
L = 2*10^(-3);
a = 0;
B = (6/(sqrt( L*C )))/(sqrt((1/(L*C))-(R/(2*L))^2));
Vc_0 = [6/(sqrt( L*C )) , 11];
odefun = @(t,Vc) [Vc(2);...
((-600)*(Vc(2)))-((50000000/3)*(Vc(1)))];
tspan = linspace(0, 0.04, 4000000);
[t,Vc] = ode45(odefun, tspan, Vc_0);
subplot(2,1,1);
plot(t,Vc)
xlabel('time(s)', 'FontSize',16, 'FontName','Arial','FontWeight','bold')
ylabel('Vc', 'FontSize',16, 'FontName','Arial','FontWeight','bold')
title('Numerical Methods', 'FontSize',16, 'FontName','Arial','FontWeight','bold')
texact=[0:0.00000001:0.04];
Vexact=exp((-300)*texact).*(a*cos(((sqrt(596760000)*texact)/6))+B*sin(((sqrt(596760000)*texact)/6)));
subplot(2,1,2);
plot(texact,Vexact)
xlabel('time(s)', 'FontSize',16, 'FontName','Arial','FontWeight','bold')
ylabel('VCexact', 'FontSize',16, 'FontName','Arial','FontWeight','bold')
title('Exact Solution', 'FontSize',16, 'FontName','Arial','FontWeight','bold')
3 个评论
What strikes me is that your exact solution does not start at (0,11). The fact that you choose alpha = 0 in the "analytical solution" automatically leads to Vc(0) = 0 which contradicts your initial value for Vc.
And the order of your initial conditions in the numerical solution is wrong:
Vc_0 = [6/(sqrt( L*C )) , 11];
I suggest that you insert the expressions depending on C, R and L for the coefficients of the ODE instead of error prone numerical values. E.g. exp((-300)*texact) must be exp(-3000*texact) instead (I didn't check the other numerical values).
Are you sure that the initial conditions for Vc and dVc/dt are given in SI units ?
Prince Nino
2024-9-20
@Prince Nino, Because you didn't use the values derived from the parameters R, L, C in the code. Instead, you entered the coefficients '600' and '50000000/3' manually.
%% Parameters
C = 3*10^(-6);
R = 12;
L = 2*10^(-3);
R/L
1/(L*C)
50000000/3
回答（2 个）
James Tursa
2024-9-19
编辑：James Tursa
2024-9-19
Assuming your state is [vc,vc'] in that order, then it appears your initial conditions are backwards and should be this instead:
Vc_0 = [11, 6/(sqrt( L*C ))];
I haven't looked any deeper into your problem yet ...
Hi @Prince Nino
Perhaps you should verify the results with your Professor. You can also change the initial condition.
%% Numerical Solution
C = 3*10^(-6);
R = 12;
L = 2*10^(-3);
f = @(t, x) [ x(2);
-(R/L)*x(2) - 1/(L*C)*x(1)];
tt = linspace(0, 0.004, 4001);
ic = [0; 6/sqrt(L*C)]; % initial condition
[t, x] = ode45(f, tt, ic);
%% Analytical Solution
a = 1;
b = R/L;
c = 1/(L*C);
rpart = b/(2*a);
ipart = sqrt(c/a - rpart^2);
alpha = ic(1);
beta = (ic(1) + ic(2))/ipart;
xt = exp(- rpart*tt).*(alpha*cos(ipart*tt) + beta*sin(ipart*tt));
%% Plot results
tl = tiledlayout(2, 1);
nexttile
plot(t, x(:,1)), grid on, title('Numerical Solution')
nexttile
plot(tt, xt), grid on, title('Analytical Solution')
xlabel(tl, 'Time')
ylabel(tl, 'x(t)')
类别
在 帮助中心 和 File Exchange 中查找有关 Ordinary Differential Equations 的更多信息
2024-9-19
2024-9-21
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
